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ABSTRACT 


An error analysis study was conducted in order to assess the currently achievable 
accuracies and the future anticipated improvements in the estimation of geopotential 
differences over intercontinental locations. Extending the ideas put forward by Colombo 
(1980), an observation/estimation scheme was proposed and studied, whereby gravity 
disturbance measurements on the Earth's surface, in caps surrounding the estimation 
points, are combined with corresponding data in caps directly over these points at the 
altitude of a low orbiting satellite, for the estimation of the geopotential difference between 
the terrestrial stations. The gravity disturbance data at altitude are inferred from GPS 
measurements made from the low orbiter to the high-altitude GPS satellites, in a multiple- 
high-single-low Satellite-to-Satellite Tracking (SST) configuration. 

The mathematical modeling required to relate the primary observables to the 
parameters to be estimated, was studied both for the terrestrial data and the data at altitude. 
Emphasis was placed on the examination of systematic effects and on the corresponding 
reductions that need to be applied to the measurements to avoid systematic errors. For the 
gravitational accelerations inferred from SST data, a mismodeling related to a centrifugal 
acceleration term was identified and corrected. Alternative formulations related to the 
sampling (or discretion) and the propagated errors arising in the truncation theory 
considerations were derived. Recurrence relations for the altitude generalized truncation 
coefficients implied by Hotine's kernel, and for the degree variances implied by a first- 
order Gauss-Markov covariance model were originally developed in this study. 

The error estimation for the geopotential differences was performed using both 
truncation theory and least-squares collocation with ring-averages, in case observations on 
the Earth's surface only are used. The error analysis indicated that with the currently 
available global geopotential model OSU89B and with gravity disturbance data in 2° caps 
surrounding the estimation points, the error of the geopotential difference arising from 
errors in the reference model and the cap data is about 23 kgal cm, for 30° station 
separation. This error is expected to reduce to about 12 kgal cm, when the lower-degree 
harmonics of the reference model are improved by the incorporation of the global GPS- 
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tracking data on Gravity Probe-B. The incorporation of gravity disturbance data at altitude 
was studied using least-squares collocation with ring-averages. It was found that for a 
low-degree (Nmax = 45) reference model, the data at the altitude of GP-B (600 km) can 
improve the geopotential difference accuracy by about 7%, as compared to the use of 
terrestrial data only. However, additional high-frequency observables at lower altitude are 
needed to achieve the results obtainable when a high-degree reference model is used, and 
to this end the gradiometer data from ARISTOTELES will provide a significant 
contribution. 
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CHAPTER I 


INTRODUCTION 


The fundamental objective of geodesy is the accurate determination of the position 
and the gravity potential of points on the surface of the Earth or in the space surrounding 
the Earth. Such information is essential to support research (as well as application) in a 
number of related disciplines such as geodynamics, geophysics and oceanography. 

Historically, geodesists have divided the problem of position and gravity potential 
determination in two parts, due to the limitations imposed both by Nature and by the 
observational techniques available. In that sense, angle and distance measurements on 
the Earth would provide through triangulation and trilateration the "horizontal" 
coordinates (<J>, X) of a point. The reference surface employed in these determinations is 
the surface of an ellipsoid of revolution. The traditional practice in such determinations 
(as opposed to modem integrated approaches) emphasizes on geometric principles, while 
gravity field information is mainly used for the reduction of the surface measurements to 
the reference ellipsoid. However, since the reference ellipsoid is a surface not physically 
realizable, the third coordinate, the height of a point with respect to this surface, had to be 
determined indirectly. 

A rather easily accessible and naturally provided surface, the Mean Sea Surface 
(MSS) or Mean Sea Level (MSL), would provide the reference surface with respect to 
which geopotential numbers and heights could be reckoned. As long as the MSL, as 
realized by tide gauge observations, was identified with a unique equipotential surface-the 
geoid, spirit leveling (a highly accurate geodetic measurement type) and gravity 
observations, would provide the geopotential number of a point. With certain 
approximations involved, geopotential numbers would yield orthometric heights, while 
gravimetry would provide the geoidal undulations (e.g. through Stokes' integral) 
required for the computation of heights with respect to the reference ellipsoid. 
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Under such operational procedures, the concepts of height and geopotential 
difference become heavily inter-related; their determination constitutes the second part of 
the problem mentioned in the beginning, where the physical properties of the Earth play a 
dominant role. The advantage of MSL is that it provides a natural connection between 
continents, enabling determination of height and geopotential differences between points 
which cannot be connected by leveling. In that sense, MSS becomes the natural 
reference surface for these determinations and establishes a world vertical datum. 
Obviously, the entire setup heavily depends on the assumption that MSS coincides with a 
unique equipotential surface of the Earth's gravity field, and on the accurate realization 
and monitoring of MSL. 

1.1 Vertical Datum Inconsistencies and the Impact of Modern Space 
Techniques 

Advances in a number of areas, that occurred during the last two decades, have 
caused geodesists to reconsider the classical procedures of vertical datum definition and 
height determination. Two main reasons are responsible for that: 

(1) It has been well recognized by now that the MSS departs from an equipotential 
surface due to the presence of the Quasi-stationary Sea Surface Topography (QSST), 
whose magnitude is on the order of a meter. The presence of QSST affects the definition 
of vertical datums in two ways: 

(a) Vertical datums established with respect to different tide gauge stations do not refer 
(in general) to the same equipotential surface, their offsets being on the order of ± V2 
kgal m (1 kgal m = 10 m 2 /s 2 ). 

(b) If in the adjustment of a vertical network, more than one tide gauge station were 
constrained to have zero elevations (e.g. the 1929 General Adjustment of the U.S. 
vertical reference system), the QSST differences at these stations will cause internal 
distortions that propagate throughout the adjusted network. 

Laskowski (1983) has studied the effects of vertical datum inconsistencies on 
various gravimetric quantities by constructing certain models for the inconsistencies. 
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using also the oceanographic estimates of the QSST derived by Lisitzin (1965). 
Laskowski (ibid) also studied the likely distortions that the overconstraining of the 1929 
adjustment of the U.S. leveling network might have caused. He concluded that the 
internal distortions are most severe near the tide gauges where incorrect zero-constraints 
were imposed. On a global basis, the gravity anomaly errors implied by the vertical 
datum inconsistency models he used, when analyzed harmonically, indicated that most of 
the power of the gravitational signatures is concentrated below spherical harmonic degree 
60. 

(2) The advent of modern space techniques, such as Very Long Baseline 
Interferometry (VLBI), Satellite Laser Ranging (SLR) and the Global Positioning System 
(GPS), has changed fundamentally the position determination procedures. 

Laser ranging to the high-altitude LAser GEOdynamic Satellite (LAGEOS) has 
enabled the determination of geocentric coordinates for a number of globally (but not 
evenly) distributed permanent tracking sites to accuracies at the ± 5 cm level (Smith et al., 
1985). VLBI measurements, on the other hand, are capable of determining baseline 
vectors between stations separated by 7000 km, accurate to a few centimeters, from 24- 
hour observing sessions (Herring, 1986). Being a geometric technique however, 
insensitive to the stations' locations with respect to the Earth's center of mass, VLBI 
requires that, at least one station's geocentric coordinates in the network, be defined from 
another source (e.g. from SLR observations). Finally, relative positioning using the 
GPS has already proved its capability to yield baseline components accurate to 2-3 ppm 
of the baseline length on a local or regional basis. With the full satellite constellation in 
orbit, dual frequency receivers and orbit determination using a fiducial network of 
tracking stations, it is expected that the system's performance will be at the 3 mm plus 
0.01 ppm of the baseline length (Carter et al., 1989). Apart of the individual accuracy 
achievements in positioning, of each of the above techniques, the recent (January 1, 
1988) establishment of the International Earth Rotation Service (IERS) (Mueller, 1988) 
enables the optimum combination of the results obtained by the above techniques (as well 
as Lunar Laser Ranging, and Doppler techniques), and the definition and maintenance of 
a Conventional Terrestrial Reference System (CTRS). The realization of such a system is 
accomplished through a network of stations whose geocentric coordinates are estimated 
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within a few centimeters and which constitute the Conventional Terrestrial Reference 
Frame (CTRF). 

The main impact of the space techniques on the problem of height determination 
comes from their ability to provide all three components of the geocentric Cartesian 
coordinate vector, and thus geometric (ellipsoidal) heights as well. Hence, one of the 
two uses of leveling (that of providing vertical position information) becomes 
unnecessary. Consequently, in the context of modem position determination techniques, 
the concepts of height and geopotential difference may be well distinguished. Substantial 
effort has been made during the last few years, by a number of researchers (Engelis et al. 
(1985), Kearsley (1986), Schwarz et al. (1987), Rapp and Kadir (1988)) to investigate 
the data requirements and the attainable accuracies in deriving orthometric height 
differences from GPS-derived ellipsoidal height differences and gravimetric undulation 
differences. Such procedures aim to eliminate the need for spirit leveling in regional 
applications, essentially reversing the traditional geodetic practice in height determination. 
The aforementioned studies indicate that orthometric height differences accurate to about 
2 ppm (of the baseline length), for lines 10-70 km long, are attainable, provided good 
gravity data coverage exists and one carefully accounts for topographic effects. Although 
such procedures may never reach the accuracy level of first-order spirit leveling, they 
certainly provide a promising cost-effective alternative for lower order vertical control on 
a regional basis. 

1.2 The Task of Vertical Datum Unification and the Proposed 

Approaches 

From the previous discussion, it becomes obvious that the internal inconsistencies 
in continental vertical networks, due to over-constraining of tide gauge stations, can be 
removed by re-adjusting the networks using only one geopotential number constraint 
(minimum constraint solution) (Laskowski, 1983). Such procedure was followed in the 
establishment of the United European Leveling Network 1973 (UELN-73), where the 
Normal Amsterdam Piel (NAP) was held fixed to an arbitrarily assigned geopotential 
number as the "origin" point of the network (Kelm, 1985). The same procedure is 
adopted for the definition and re-adjustment of the North American Vertical Datum 1988 
(NAVD 88) (Zilkoski, 1986). 
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The removal of internal distortions from continental vertical networks does not 
pose any theoretical problems, nor does it require any additional observations other than 
the geopotential differences obtained from leveling and gravimetry. The implications 
however (e.g. in map-making) of changing the height system for an entire country or 
continent have to be considered carefully. In addition, even if a tide gauge were to be 
selected as the "origin" point of such network, the presence of the QSST would cause 
some points along the coast to be below, and some above, the MSL (Colombo, 1985a). 
Land "below sea level", shown on a coastal map, would be undesirable from the practical 
point of view; however, as Colombo (ibid) pointed out, such problems may be avoided 
using a local vertical reference for the region, based on a local tide gauge. Provided that 
at least one station is connected to both the local and the national (or international) 
network, both sets of "heights" may be converted to each other unambiguously. Such 
procedures are conceptually very similar to the use of local "best-fitting" ellipsoids as 
opposed to a global "Mean-Earth" ellipsoid (Heiskanen and Moritz, 1967, section 5-11) 
in "horizontal" network applications. 

According to the above, the definition and realization of a global vertical network 
finally reduces to the task of accurate determination of geopotential differences between 
points that cannot be connected by spirit leveling (combined with gravity observations), 
i.e. between points separated by ocean. Once these geopotential differences are 
determined, e.g. between the "origins" of the various (internally consistent) continental 
networks, an arbitrary value for the potential at one "origin" point, is enough to provide a 
global uniform system of geopotential numbers. 

Although the need for internally consistent continental networks arises from 
practical considerations, establishing intercontinental vertical connections is primarily of 
scientific, rather than operational interest, and the necessity of such connections has been 
debated in the literature. One advantage that a global vertical network possesses, is that it 
will enable referencing of regional gravity anomaly databases to a uniform system, thus 
freeing global anomaly fields from regional systematic errors induced by height 
inconsistencies (Rapp, 1983; Heck, 1990). Colombo (1985b) however, has pointed out 
that GPS observations and gravimetry may enable one to switch from gravity anomalies 
to gravity disturbances, in future gravity surveys, the latter ones being independent of 
orthometric height inconsistencies. Although this is true, it is also unreasonable to 
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believe that the wealth of gravity anomaly data acquired in the past will be simply 
abandoned under such operational procedures. 

The approaches proposed for the unification of vertical datums are based on either 
one of the two following principles: 

(a) Since the MSS does not coincide with an equipotential, define and realize in some 
manner another equipotential surface to be the reference; if both the MSS as well as its 
departures from the reference surface are accurately determined and monitored, the 
transoceanic connections based on tide gauges can be maintained. 

(b) Abandon MSS as a transoceanic connection and seek alternative techniques of 
estimating accurately geopotential differences between points separated by ocean. 

Colombo (1985b) classified the observational techniques and estimation 
procedures proposed up to now for the unification of vertical datums into four main 
categories; the first two follow principle (a) above, while the other two (b). 

(1) Oceanog raphic Approach 

This technique, outlined by Cartwright (1985) considers the reference surface to 
be at a fixed geopotential number above an isobaric (i.e. equipressure) surface, 2000 
decibars below the annual MSS at some specified epoch. Cartwright's technique 
employes hydrographic measurements along profiles that extend from shallow areas near 
a tide gauge, to ocean depths of more than 2 km. These profiles are selected to coincide 
with the repeat groundtracks of an altimeter satellite that can provide estimates of the 
ellipsoidal height of the sea surface. Steric and geostrophic leveling are used to determine 
the shallow point elevation relative to the deep isobaric surface, while pressure 
measurements and spirit leveling are used to connect the shallow point to the nearby tide 
gauge. 

From the theoretical point of view, arguments against this technique have been 
raised, related to the accuracy and suitability of the isobaric surface ("level of no motion") 
hypothesis (Colombo, 1985b). In addition, the geostrophic leveling used to connect the 
deep ocean location to the shallow one, may not be accurate enough to model the complex 
ocean dynamics near the continental boundary (Wunsch and Gaposchkm, 1980). From 
the practical point of view, several oceanographic restrictions limit the selection of sites 
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where the technique could be applied (e.g. sites must have narrow continental shelves, 
should not be affected by river discharge, strong currents etc.). In addition, at least a full 
year of (nearly) simultaneous observations at all sites is required to establish the 
transoceanic links, and the resulting vertical datum would have a rather strong time 
dependency and would require periodic re-definition. Although periodic maintenance 
would be necessary for a global vertical datum, no matter what technique is used to 
define and realize it (since the Earth's gravity field undergoes secular and periodic 
changes due to a variety of geodynamic phenomena such as post-glacial rebound, mass 
redistribution etc.), alternative techniques may offer better temporal stability. 

(2) Altimetrv-Gravimetrv Approach 

This technique, considered variously by Mather et al. (1978) and Wunsch and 
Gaposchkin (1980), utilizes altimetric observations in combination with ocean gravimetry 
and aims to the simultaneous recovery of the QSST and the geoidal undulation. The 
estimation technique proposed by Wunsch and Gaposchkin (ibid) is in essence least- 
squares collocation, and the separation of the geoidal from the QSST signal is aided by 
the use of prior information in the form of a-priori degree variances for these signals. 
One of the limitations of such procedures comes from the inaccurate and insufficiently 
sampled oceanic gravity measurements. However, the ideas of Wunsch and Gaposchkin 
(ibid) have been pursued further by a number of investigators (Wagner, 1986; Engelis, 
1987) and recent analyses (e.g. Denker and Rapp, 1990) have verified the ability of 
satellite altimetry to determine the long-wavelength features of the global ocean circulation 
and simultaneously provide improved estimates of the long-wavelength part of the 
oceanic geoid, as well as, improved orbital parameters. Such solutions though, aim to 
resolve global features of the QSST, with wavelengths greater than about 4000 km; 
detailed local determinations in the shallow areas near the tide gauge locations are further 
complicated by the inaccuracies of the tidal models. 


(3) Geodetic Approach 

The main ideas in this direction have been put forward by Colombo (1980). He 
considered the realization of a global vertical network by a set of benchmarks whose 
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geocentric coordinates and geopotential differences are accurately determined. In his 
analysis, SLR techniques would provide the geocentric coordinates, while surface 
gravimetry, spirit leveling and a low-degree geopotential model would provide the 
geopotential differences, in a least-squares collocation solution. He has also considered a 
modification of the observed gravity anomalies to avoid possible contamination of the 
estimated geopotential differences, by systematic errors in the anomalies due to height 
inconsistencies. Additional error analysis performed by Hajela (1983) has indicated that 
Colombo’s technique could provide a US - Europe connection at that time, to an accuracy 
of about 50 to 60 kgalcm. The reference geopotential model used at that time was the 
OSU81 (Rapp, 1981a). The disadvantage of Colombo’s technique lies on its dependence 
on the reference geopotential model. That is, the transoceanic “connection” of stations is 
carried out analytically, through the long-wavelength information provided by the 
reference model. Consequently, any errors of the low-degree harmonics propagate 
directly to the estimated geopotential differences. This has also been manifested in the 
contribution of the errors of the model to the total error budget, in the error analyses 
performed. 

Hein and Eissfeller (1985) have discussed the application of an integrated geodesy 
adjustment approach to the problem, considering stations at (or near) tide gauges and 
additional observations to those considered by Colombo (1980), such as altimeter 
measurements and deflections of the vertical. Brovar (1988) and Rummel and Teunissen 
(1988), on the other hand, have considered the analytical modification to the solution of 
the geodetic boundary value problem, so that vertical datum offsets can be introduced as 
unknowns and estimated by solving a linear system of equations. The technique of 
Rummel and Teunissen (ibid) requires the same type of observables as the one of 
Colombo (1980). However, the former, is in the strict sense applicable only if vertical 
datum offsets are identified worldwide and introduced as unknowns in the linear system 
to be solved, provided also that global coverage of gravity anomalies is available; such 
restrictions have been circumvented in Colombo’s (ibid) approach. 

(4) Relativistic Approach 

Based on a different physical principle, this technique, introduced by Bjerhammar 
(1985), aims to the direct measurement of geopotential differences, using the effect of 
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gravitation on frequency standards, as predicted by the theory of general relativity. 
Unfortunately, this effect is too small to be detectable by the frequency standards 
available at present, and application of this technique has to be postponed until orders of 
magnitude improvements are achieved in the accuracy of frequency standards. 

The scientific and technological advances that have been achieved, or are 
anticipated for the near future, appear to favor at present the implementation of 
oceanographic or geodetic type of approaches for vertical datum connections. Both 
techniques (but more critically the oceanographic) require some internationally co- 
ordinated observation campaign for their implementation (Rapp, 1987); given the limited 
interest (from the operational point of view) for unification of vertical datums, techniques 
that promise best results while taking maximum advantage of data that either exist, or will 
be collected for other investigations, should be preferred. Also, preference should be 
given to those techniques that provide longer temporal stability in the resulting unified 
global vertical network, and thus require less often re-definition and maintenance. 

1.3 Motivation, Objective and Organization of the Present Study 

Two types of geodetic projects currently under planning and/or development create 
a favorable situation for the implementation of geodetic techniques for vertical datum 
connections. 

First, the incorporation of GPS receivers on-board a number of satellites 
scheduled for launch in the 1992-1995 time frame. Table 1 (see also (Colombo, 1990)) 
provides some information related to these missions. With the GPS receiver on-board 
the lower orbiter observing simultaneously as many as seven GPS satellites, a Multiple- 
High-Single-Low Satellite-to-Satellite Tracking (SST) configuration is established (Jekeli 
and Upadhyay, 1990). Such configuration enables one to estimate all three components 
of the total inertial acceleration at the altitude of the low orbiter. Provided non- 
gravitational accelerations can either be measured (as is proposed for ARISTOTELES), 
or effectively be eliminated (as in the case of the drug-free GP-B spacecraft), the result of 
such an observational system may be viewed as a dense grid of uniform quality 
"observations" of the gravitational acceleration vector at the altitude of the low orbiter. 
There are two ways that such "observations" can contribute to the solution of the problem 
at hand: 
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(a) Analysis of the global set of measurements collected during the lifetime of the low 
orbiters, is capable of producing global geopotential models whose quality is orders of 
magnitude better than current state-of-the-art models, and whose resolution (half 
wavelength) ranges from about 600 km (TOPEX/POSEIDON) to 100 km 
(ARISTOTELES). These figures have been assessed initially for GP-B by Smith et al. 
(1988) using an analytical approach, and have been verified recently through more 
elaborate and complete simulation studies (Pavlis, E. et al., 1990). Given the 
importance of an accurate geopotential model in the implementation of geodetic 
approaches for vertical datum connections, a major contribution to the achievable 
accuracies is to be expected from the above missions. 

(b) The dense grid of observations at altitude may help resolve localized signatures of the 
gravitational field on the surface of the Earth, provided the satellite is low enough and/or 
instrumentation on-board provides additional measurements, more sensitive to the finer 
details of the field (e.g. the gravity gradiometer in the case of ARISTOTELES). 


Table 1. Future Satellite Missions Expected to Carry GPS Receivers On-Board. 


Name 

Scheduled 
Launch Date 

Altitude/ 

Inclination 

Description 

TOPEX/ 

POSEIDON 

1992 

-1336 km 
66.02° 

USA and French Altimeters and other 
oceanographic instruments. GPS 
receiver under development. 

GRAVITY 

PROBE-B 

(GP-B) 

1995 

~ 600 km 
90° 

Stanford General Relativity Gyroscop- 
ic Experiment Drag-Free. 2-year 
mission. 

ARISTOTELES 

1995 

-200 km 
96.3° 

European Space Agency's Gravity 
Gradiometer. 6-month mission. 
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A second type of geodetic activity, pertinent to the vertical datum connection 
problem, is the geodetic fixing of Tide Gauge Bench Marks (TGBMs) (Carter et al., 
1989), a project that has been initiated by the International Association for Physical 
Sciences of the Ocean (IAPSO). Highly accurate determination of the geocentric 
positions of the TGBMs (as well as accurate monitoring of their motions) is necessary, in 
order to separate the crustal motions of the TGBMs from the apparent changes of the 
MSL observed at their locations, and thus, provide the means of investigating the 
possibility that global MSL shows a rising trend resulting from global warming (IAPSO, 
1985). Although monitoring changes in the global sea level does require 
establishment of a global vertical network, the data that are to be collected for the former 
purpose, can be of use for the latter as well. In addition, if geopotential differences 
between TGBMs can be determined accurately, without any reliance to MSL, they can be 
used to provide independent control over oceanographically determined QSST differences 
between the TGBMs. 

The above developments provided the motivation to reconsider (from the geodetic 
point of view) the problem of estimating geopotential differences between points 
separated by ocean. As mentioned already, this is equivalent to the problem of 
connecting different vertical datums to a common reference, thus establishing a uniform 
global vertical datum. Following the ideas of Colombo (1980), a "global vertical datum 
is defined here as a network of benchmarks situated on various continents, between 
which a set of estimated geopotential differences is given. The technique used to estimate 
these geopotential differences is based on an observational system that attempts to 
improve the one considered by Colombo (ibid), and overcome some of its limitations. 
The basic "components" of the current observational setup are shown in Figure 1. 

Two benchmarks, BMA and BMB, are considered to be situated on the same or on 
different continents (separated by ocean). These are VLBI or SLR stations, whose 
geocentric coordinates are known to sub-decimeter accuracy, and in addition they are 
equipped with calibrated absolute gravimeters, so that the magnitude of absolute gravity 
is also known accurately at these sites. Differential GPS observations, and relative 
gravity measurements, provide estimates of the gravity disturbance at points inside "caps" 



Figure 1. The Basic Observational Geometry. 

centered at the corresponding benchmarks. Consider also, "caps" centered directly over 
these benchmarks, at the altitude of a satellite which carries a GPS receiver on-board. 
Inside these caps, at satellite locations whose geocentric coordinates are known from the 
GPS tracking, some functional of the gravitational potential has been estimated (e.g. the 
three components of the gravitational acceleration vector). 

From the data requirement point of view, apart of information that either exists 
already or will be obtained as part of other geodetic activities, the additional observations 
required here are the absolute gravity measurements at the benchmarks and the gravity 
disturbances in the caps surrounding them. Careful coordination (in terms of site 
selection) with the project of geodetic fixing of TGBMs (Carter et al., 1989) can reduce 
the number of additional absolute gravity measurements required, while kinematic (or 
rapid static) GPS techniques can provide an efficient way of performing the gravimetric 
surveys inside the caps. 

Based on the above information, using also a reference geopotential model, one 
can estimate the potential difference between the benchmarks BMA and BMB, using 
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least-squares collocation (Moritz, 1980). In addition, such estimation technique, enables 
one to derive measures of the accuracy of the resulting estimates. The main objective of 
this study is to provide accuracy estimates for the resulting geopotential differences, 
based on realistic assumptions for the errors associated with the input data, and 
considering the observational setup of Figure 1, or variations of it. More precisely, the 
error analysis to be presented, will address the following issues: 

(1) The attainable accuracies for the geopotential differences, if state-of-the-art 
geopotential models, developed in the absence of the anticipated missions, are to be 
used in the estimation. 

(2) The corresponding accuracies using models of the quality expected to be obtained 
from various future missions. 

(3) The contribution of observations at altitude to the estimation of the geopotential 
differences, and the possible improvements in accuracy through the incorporation 
of such observations. 

It is worth noticing that the procedure considered here does not involve spirit 
leveling at all. In more realistic configurations, involving more than one benchmark per 
continent, if geopotential differences between these benchmarks are available from 
leveling (combined with gravimetry), they could (and should) be used in a simultaneous 
adjustment, to increase the redundancy and strengthen the solution (as in the estimation 
scheme described by Colombo (1980)). In that sense the current error analysis 
corresponds to a worst-case scenario. 



CHAPTER H 


MODELING ASPECTS 

The accuracy that can be achieved in the estimation of the geopotential differences 
using the configuration described in section (1.3), depends on the accuracies of the 
primary observables involved in the estimation, and on the way according to which 
random and systematic errors in the observed quantities propagate to the estimated 
values. The former can be assessed from information related to the performance of the 
sensing instruments involved in the data acquisition process (e.g. gravimeters); the latter 
requires the development of analytical formulations that relate the primary observables to 
the quantities of interest, and constitutes the subject of this chapter. More precisely, the 
following paragraphs focus on the modeling of the gravity disturbance observations on 
the surface of the Earth, and of the components of the gravitational acceleration vector at 
altitude as obtained from a Satellite-to-Satellite tracking configuration. 

2.1 The Boundary Condition Implied by Gravity Disturbance 

Observations 

Let (r, 0, X) denote geocentric distance, geocentric co-latitude and longitude 
respectively. The following notation definitions are adopted: 

V(r, 0, X) : true gravitational potential of the Earth. 

V e (r, 0) : gravitational potential of a reference ellipsoid of revolution whose 
surface is an equipotential surface of its gravity field. 

V m (r, 0, X) : true value of the gravitational potential of the Earth, that arises from 
all the harmonics only up to maximum degree M. 

d>(r, 0) : true centrifugal potential of the Earth. 


14 


It is well known (Heiskanen and Moritz, 1967, p. 47) that: 
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<p(r, e) = ^o 2 r 2 sin 2 0 (2J) 

where 0) is the spin rate of the Earth. An estimate of <X>(r, 0) may be obtained from 
estimated values of r, 0 and co. Assuming for a moment that (r, 0) are perfectly known, 
the error in the estimated value of 0(r, 0) due to an error e® in © will be: 

e8(r, o)© = cor 2 sin 2 0e(o (2.2) 

Using the nominal values r = 6371 km, co = 7292115 x 10 _11 rad/s, and e® = 0.1 x 10*H 

rad/s (Chovitz, 1988), one has for a point on the equator (where £<1>® becomes 
maximum): 


maxEO (r, 0 ) 0 , » 3 x 10' 3 m 2 s' 2 = 3 x 10^ kgal m 


(2.3) 


Such error in potential translates to linear ("height") error of about 0.3 mm and therefore 
is, in the context of this study, negligible. Hence, in all subsequent analysis, the spin 
rate of the actual Earth will be considered perfectly known and equal to the spin rate value 
used in the definition of the reference ellipsoid, i.e. 

e$ (r, 0)o = 0 , (2.4) 

and the centrifugal potential estimated from perfectly known (r, 0) but approximately 
known CD, will be considered identical to the true centrifugal potential of the actual Earth. 
The following notation will be used: 

W(r, 0, X) = V(r, 0, X) + <D(r, 0) (2 .5) 

U e (r, 0) = V e (r, 0) + $>(r, 0) (2 .6) 


U m (r, 0, X) = V m (r, 0, X) + <D(r, 0) 


(2.7) 
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so that, W represents the true gravity potential of the Earth, U e the gravity potential of the 
reference ellipsoid, and U m the value of the part of the gravity potential of the Earth 
that contains harmonics of the true gravitational potential only up to maximum degree M. 

*n\ in , , 

V m and U m should not be confused with estimates V and U of these quantities, 
obtained e.g. through satellite perturbation analysis; the former are true values while the 
latter are contaminated by the commission error of the estimated harmonic coefficients. 

The disturbing potential T at a point P(r, 0, X), is now defined by Heiskanen and Moritz 


(1967, section 2-13): 

T P = W P - (2.8) 

with respect to the ellipsoidal field, and by: 

Tp = W P - Up (2.9) 

with respect to the truncated field U m . Due to the previous definitions and the 
assumption made concerning the centrifugal potential, one has: 

T P = Vp - Vp (2.10) 

Tp? = V P - Vp 1 (2.11) 

In addition, the following quantities are introduced (Heiskanen and Moritz, ibid): 

gravity vector : g p = grad(W)p (2.12) 

normal gravity vector : yp = grad(U e )p (2. 1 3) 

gravity disturbance vector : 5gp = gp-yp (2.14) 


and the geometrical relations associated with these definitions are illustrated in Figure 2. 
According to the previous definitions and due to linearity of the gradient operator one has: 


5gp = grad(T)p 


(2.15) 
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e h 



Figure 2. Geometry Associated With the Gravity Disturbance. 

Now, proceeding along the same lines as Moritz (1980), if x is the arc-length along the 
isozenithal (ibid, p. 345), one can derive from equation (2.15): 



(2.16) 


f — > 

where ogp is the component of 8gp in the downward direction of the isozenithal passing 
through P. The very small curvature of the normal plumb line (Moritz, 1983, p. 7) 
justifies the approximation that the normal plumb line coincides with the straight 
ellipsoidal normal, in which case the isozenithal will also coincide with the straight 
ellipsoidal normal (Moritz, 1980, pp. 345-346). Under such an assumption, equation 

(2.16) becomes: 



(2.17) 
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and 5gp will be hereon understood to represent the component of 8gp along the 
downward direction of the straight ellipsoidal normal through P. If e^ is the unit vector 
in the direction of increasing ellipsoidal height, and e n the unit vector along the normal 

plumb line (pointing outwards), then the assumption made before implies: 



(2.18) 


and equation (2. 17) may be written as: 


5g p = (gP-7p)*(-e h ) (2.19) 

where the dot denotes scalar product. According to (2.18) the last equation becomes: 


Sgp = * gp • e h - lypl 

= - grad(W)p • e h - lypl or finally: 


- 1 Ypl = 8g p + 



( 2 . 20 ) 


On the other hand, if e H is the unit vector along the direction of increasing orthometric 
height, one has: 


Igpl = grad(W) P ■ (- e H ) 


or, (Pavlis, 1988, equation 2.31) : 






aw 

T1 

M3tp 


aw 

NcostpaX. /p 


( 2 . 21 ) 


( 2 . 22 ) 
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where, 0 is the total deflection of the vertical, £ and T| its latitudinal and longitudinal 
components respectively (Heiskanen and Moritz, 1967, p. 83), M and N the meridional 
and prime vertical radii of curvature respectively and <p denotes geodetic latitude. 

Adding equations (2.20) and (2.22) by parts, one obtains: 

Igpl - ly P l = 5g p + e p (2.23) 


where: 


e P = 


i^ ft dW 

(1 - cos 0/—- £ 




aw 


8h M9<p Ncos cp3A.Jp 


(2.24) 


From equations (2.17) and (2.23) one has: 


Igpl - 1 y P l = - 


aT| 

i3hj 


+e P 


(2.25) 


However, it was shown by Pavlis (1988, equation 2.52), that: 


r a°] 

,3h >p 


— | - e 2 sin0pcos0p(-^— ) +0(e 4 ) 
\dr/p \ra0/ P 


(2.26) 


where e is the first eccentricity of the reference ellipsoid. Hence, neglecting terms of the 
order e 4 and higher, one has from (2.25): 


Igpl - |y p | = - |— -j +(£h)p+Ep 


(2.27) 


where. 


( £ h)p = 


e 2 sin0cos0 


3T 


\r90/Jp 


(2.28) 
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Equations (2.25) or (2.27) represent the boundary condition of a Neumann-type 
boundary value problem (bvp) of potential theory (Kellogg, 1954), linearized with 
respect to a Somigliana-Pizzetti normal gravity field (Heiskanen and Moritz, 1967, 
section 2-7). From equations (2.22) and (2.24) follows that: 


lgpl = - 


aw\ 


\ dh )p 


+ Ep 


(2.29) 


Equation (2.25), when compared with (2.29), indicates that linearization with respect to 
the aforementioned normal field, removes the centrifugal terms from (2.29). In addition, 
linearization permits the truncation of the Taylor series expansion of 3*/3h, around d'fdr, 
to terms of 0(e 2 ) in (2.27). Provided that the effect of the mass of the atmosphere has 
been taken into account analytically (see section 2.3), the disturbing potential T may be 
considered harmonic outside the topographic surface (S). The problem at hand then, is to 
determine the function T, harmonic outside (S), and regular at infinity, subject to the 
boundary condition (2.25) or (2.27), both valid over the known boundary (S). This is a 
fixed bvp; however, both (2.25) and (2.27) represent oblique derivative boundary 
conditions, since neither the ellipsoidal normal, nor the geocentric radius vector, are 
normal to the surface (S), where the boundary values are given. Equation (2.25) 
contains the effects of the approximation (2.18), while (2.27) contains in addition the 
effects of the neglected terms of 0(e 4 ) and higher in the Taylor series expansion of 979h 
around d'/dr. Equation (2.27) may also be written as: 


lgpl-ly p l-[(£h)p+ep]=- 


r 3T 


(2.30) 


or, due to equations (2.10) and (2.1 1) : 

Ta(v">-V e ) 


Igpl - lypl - [(£h)p + £p] + 


dr 


far”! 

dr 


(2.31) 


For the purpose of future reference, equation (2.25) is also repeated here, written in the 
form: 


Igpl - lypl - e P = - 



(2.32) 
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Regardless of the form of the boundary condition that one adopts, the primary 
observable from surface gravimetry is the magnitude of the gravity acceleration at the 
surface point P. The magnitude of the normal gravity at the same point can be computed 
to any degree of precision, once the reference ellipsoidal gravity field is defined and the 
geocentric coordinates of P are given (Heiskanen and Moritz, 1967, section 6-2). Hence, 
the gravity disturbance 

5g P = Igpl - lypl , (2.33) 

as obtained from surface gravimetry, contains the observation errors of Igpl and the errors 
in lypl induced by the errors in the geocentric coordinates of P. The latter is a 
misregistration error, i.e., the actual observation refers to a different location than the one 
defined by the geocentric positioning results. 

From the point of view of an analytical formulation for the solution of the current 
bvp, it is obvious that 5g is related to the unknown disturbing potential, in a rather 
complicated manner, due to the presence of the atmospheric effects, the ellipsoidal terms 
(Ej|)p and ep, and even more critically due to the fact that it is defined over the 
topography, a very complex surface which cannot be describe analytically (Holota, 
1985). Fortunately however, the dominant spherical character of both the shape, as well 
as the gravity field of the Earth, and the fact that the atmospheric mass amounts to only 
about 10' 6 (Moritz and Mueller, 1987, p. 4) of the mass of the Earth (hence the 
atmospheric attraction is about 10" 6 the attraction of the rest of the Earth's mass), have the 
consequence that the solution of the current bvp can be approximated to a high degree of 
accuracy by the solution of a second bvp of potential theory for the space exterior to a 
sphere. The latter solution can subsequently be refined by appropriate corrections that 
account for the differences between the real world and the idealized situation. The 
analytical solution of Neumann's bvp for the exterior of a sphere is the subject of the next 
section, while the refinements to this solution are considered afterwards. 

2.2 Solution of Neumann's Boundary Value Problem for the Exterior 

Space of a Sphere 

In the general case, the statement of the exterior Neumann's problem is given as 
(Kellogg, 1954, p. 246): 
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"Determine a function U, harmonic in the infinite region outside a closed surface 
S, if its normal derivatives 3U/3n assume on the surface S prescribed values". 
Harmonicity over an infinite region will be understood to include the demand for 
regularity at infinity (ibid, p. 217). Such demand ensures uniqueness for the solution of 
the exterior problem, unlike the interior one, where the unknown function can only be 
determined up to an additive constant (ibid, p. 246). 

To develop an integral formula for the solution of the exterior problem, one must 
find a Green's function of the second kind G(P ; Q), such that: 

iG(p ;Q )^«3ds Q 

S (2.34) 

where, P is outside S and Q defines the location of the variable area-element dS Q on S. 
In the current case the surface S will be the surface of a sphere of radius R, centered at 
the origin O of the coordinate system used (see Figure 3). 



Figure 3. Geometric Relations Used in the Derivation of Hotine's Kernel. 
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The unknown function U will be denoted here T, while the boundary values 9T/9r with 
opposite sig n will be denoted 5g, i.e. 

8g = -^ (2.35) 


on the sphere (O, R). For the determination of the required Green's function one may 
proceed as in Hotine (1969) (see also (Sjoberg, 1990)). From Figure 3 one has: 


£2 = r 2 + r2 . 2Rt cosy 


(2.36) 


and setting: 



(2.37) 


yields: 


L = (l + k 2 - 2k cosy)' 1 / 2 

l v ; (2.38) 

The right-hand side of the last equation is readily recognized to be the generating function 
of the Legendre polynomials (Davis, 1975, p. 365); hence: 


f-= X k n Pn(cOS \|/) 
« n=0 


(2.39) 


where, P n (t) denotes the n^-degree Legendre polynomial of the first kind with argument 
t (Heiskanen and Moritz, 1967, section 1-11), and the infinite series in (2.39) is 
absolutely and uniformly convergent for k<l. Considering k as an independent variable 
in (2.39) and integrating both sides of (2.39) with respect to k, one obtains: 



P n (cos y) . 


(2.40) 



Now, as Hotine (1969, p. 31 1) has observed: 
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2n + 1 - ? - 1 

n+1 n+1 


so that, from (2.39) and (2.40), taking also into account (2.37), one has: 



dk = ^ k n+1 P n (cos \|/) . 

* n=0 


(2.41) 


(2.42) 


The indefinite integral in (2.42) may be evaluated in closed form as: 


I 


J-dk = ln2 + ln(£ + k - cos 



(2.43) 


Considering the limiting value of the expression on the right-hand side of equation (2.43) 
as r tends to infinity, it may be easily verified that: 



dk = lnl 


(L + k- 


COS \| f 


1 - cosy 


(2.44) 


In addition. 


lim | 0 , 

r->°°U | 


(2.45) 


so that equation (2.42) implies: 


2 ^ - lnl 


£ + k - cos \jr 


1 - cosy ) 



(2.46) 


Since, 


|P n (t)|£ 1 ; -1 £t£ 1 , n = 0, 1,2,... 


(2.47) 
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one has: 


lim ( V 2n + 1 k n+1 P n (cos \j/) 

r— I n n + 1 
\n=0 


-2 

n=0 


P n (cos y) • lim {k n+1 } = 0 , 
n + 1 r->«* ' ' 


as long as k < 1. Hence, equation (2.46) finally becomes: 


(2.48) 


H(k, y) = X 2n T L j I k n+1 P n(cos y) 

O n i a 

where: 

and: D 2 (k, y) = 1 - 2kcosy + k 2 


(2.49) 


(2.50) 

(2.51) 


The series in (2.49) has the same convergence properties as the one in (2.39), for k < 1. 
The case k — 1 requires special consideration in order to establish conditional 
convergence, except of course for y = 0 (Hotine, ibid). The function H(k, y) will be 
designated hereon Hotine's function (or kernel); more precisely the form (2.50) will be 
referred to as Pizzetti's extension of Hotine's kernel, while the term Hotine's kernel will 
be reserved for the special case when k = 1 . 


The unknown function T is now identified to represent the disturbing potential 
(equations (2.8) and (2.10)). Assuming that T is harmonic outside the sphere (O, R), the 
following relation will be adopted for its expression in terms of solid spherical 
harmonics: 


T<r,e,X) = ^X(f-f X CnmYnmM) 

n=0 m— n 


(2.52) 


where: 


GM is the geocentric gravitational constant 
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Cnm arc fully-normalized, unitless, disturbing potential coefficients (with 
respect to the ellipsoidal gravitational field) 

and: 


Ynm(® » X) — P njm|(cOS 0) 


cos mX if m !> 0 
sin |rr|X if m<0 


(2.53) 


where P n m(t) is the fully-normalized associated Legendre function of the first kind with 
argument t (Heiskanen and Moritz, 1967, section 1-1 1). Denoting the n^ 1 -degree surface 
spherical harmonic (ibid, section 1-10) of T by T n (0, X), i.e., 


l(r ) 0 ) X) = GMf; CmYJq, X) 

n=0 m— n 


(2.54) 


equation (2.52) becomes: 


i(r.8,x)= S (f.jr'TJte, X) 

n-0 r (2.55) 

and both (2.52) and (2.55) are convergent for r > R (ibid, section 1-16). According to 
the last equation one has: 

& n=o (2.56) 

and thus, according to (2.35): 


8g(r. 8. X) = f X (n + X) ■ 


n=0 


(2.57) 
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Note that 8 g here is used merely as a symbol for - 3T/8r, with no explicit connection to 
the quantities that are actually observed, in contrast to the discussions of section 2. 1 . The 
relation between the current 5g and the observables will be considered later. 

Assuming that 8 g(r = R, 0, A,), on the surface of the sphere (O, R), can be expressed as a 
convergent series of surface spherical harmonics, 5g n (0, A), one has: 


8 g(R, 0,A)= £ 5g n (e,A) • 

11=0 


(2.58) 


and, due to (2.57) and (2.58), it follows that: 

Sgj(e,».)- « L | J -T„(e.x) . 


(2.59) 


The surface spherical harmonics 8g„(0, A), may be determined by (Heiskanen and 
Moritz, equation 1-71) : 


8g n (0, A)= 2iL±l|| 8g{R,e , ,A')p n (cos V )da , 
a 

where a is the surface of the unit sphere, and : 
cos\|r = cos 0 cos 0 ' + sin0sin0'cos (A - A') 

da = sin 0 'd 0 'dA' 

Equation (2.60), due to (2.59), yields : 

Tn ( 0 , A) = 2il ±_ljj 5 gj R , 0 ', A')p n (cos => (due to (2.55)) 


(2.60) 


(2.61) 


(2.62) 
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X gT*# p "( cos v) 


L n=0 


8g(R,e'A')da , 


(2.63) 


where in the last equation the orders of integration and summation have been 
interchanged. Due to (2.49), and considering (2.37), the last equation finally becomes: 


i<r, e, a.)= JLjj h(k. , vjSgfR, e', v)a< 


(2.64) 


The last equation is the desired integral formula, that determines the value of the 
harmonic function T, at any point P(r, 0, X) outside the sphere (O, R), from the values of 
its normal (radial) derivative, given continuously over the surface of this sphere. In the 
limiting case where r— »R, it can be shown easily that: 


H(l, y) ■ H(\j/) = esc y - ln(l + CSC y) 


(2.65) 


and : H(y) = £ 2 a ± 1 p n ( cos y) , 

n=0 


( 2 . 66 ) 


while equation (2.64) becomes : 


^=%lf H(v)Sg(e', r)do , 


where the constant radius R was omitted from the notation of the kernel H(y). 


From equations (2.55) and (2.59) one also obtains: 


(2.67) 


l(r, 0, *) = r£ -| T (fp 1 8g,,(e. *•) 

n=0 1 


( 2 . 68 ) 
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which represents the solution of the bvp in question, in terms of spherical harmonics 
(compare with (Heiskanen and Moritz, 1967, equation 1-91)); the surface spherical 
harmonics 5g n (0, X) are again obtained from equation (2.60). 

Comparing equation (2.64) to (2.34), taking also into account the definition (2.35), it 
becomes obvious that the Green's function of the second kind G(P ; Q) for the problem at 
hand is : 


G(P;Q) = -Rh(& y) 


(2.69) 


The function H(k, \j /) possesses a number of properties which are given next : 

(1) H(k, y)>0 ; 0<k<l, Ocy^Tt (2.70) 

This can be proved easily if one recognizes that for the above range of k and \|/, it holds 
true that : 


2k ^ D + k - 1 
D a 1 - t 


(2.71) 


and utilizes a series expression for the quantity ln(l + x) for x > 1/2. 

(2) V 2 H(k, y) = 0 ; k < 1 , (2.72) 

i.e., H(k, y) is harmonic outside the sphere (O, R), as it can be seen immediately from 
the expression (2.49). 

(3) - r S^U*(^ R2 > , (2.73) 

3r t 

which may be proved either by direct differentiation performed in equation (2.50), or, 
more simply, by differentiation of the series expression (2.49) utilizing also the relation 
(Heiskanen and Moritz, 1967, p. 35) : 

r n=0 


(2.74) 
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Note that the right-hand side of (2.73) is exactly the kernel function of Poisson's integral 
(ibid, equation 1-89). In addition, from (2.73) may be easily verified that: 

9H(k, v) \ _ 0 . 

31 I (2.75) 

(the case = 0 requires use of L'Hospital's rules to prove the last relation). This relation 
is actually a boundary condition that is imposed on the Green's function of the second 
kind, developed for the solution of the second bvp of potential theory (Roach, 1970, 
equation 9.88), specialized here for the case of spherical boundary. 

A graph of the function H(\jr) (equation (2.65)), is given in Figure 4; for 
comparison purposes, the figure also illustrates Stokes' function S(\j/) (Heiskanen and 
Moritz, 1967, equation 2-164). 

Although gravity disturbances are geometrically and conceptually simpler than 
gravity anomalies, traditional geodetic practice has relied heavily on the latter, for the 
determination of the external potential of the Earth. The underlying reason, is the 
requirement for a known boundary for the definition of the disturbances, unlike the 
anomalies (Hotine, 1969, p. 314). Using observations of the gravity anomaly vector and 
the gravity potential, all over the unknown surface of the Earth, Molodensky's bvp is 
formulated as a non-linear, free bvp, whose solution determines not only the external 
potential, but also the physical surface of the Earth. Linearization in that case requires the 
introduction of both a normal gravity field and an auxiliary known surface (the telluroid), 
in order to reduce the original problem to a linear, fixed bvp, and thus enable a tractable 
solution. The definition of the telluroid requires additional conditions to be imposed on 
either the potential anomaly AW (and the directions of actual and normal plumb lines) 
(Marussi mapping) or the gravity anomaly vector Ag (gravimetric mapping) (ibid, pp. 
338-339). In the current case, the physical surface of the Earth is considered known, 
hence the introduction of the telluroid (and subsequently of the above conditions on AW 
or Ag), is unnecessary. The normal gravity field, is introduced in order to remove 
centrifugal terms from equation (2.29) and enable an early truncation of series related to 
ellipsoidal terms. 


lim 

r->R 
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Finally, from the previous derivations, it is obvious that the use of gravity 
disturbances circumvents the problems related to the inadmissible harmonics of degrees 
n = (0, 1), which have to be suppressed from the disturbing potential for the derivation of 
Stokes' integral (Hotine, 1969, pp. 317-318). While global gravity anomaly data, input 
to Stokes’ integral should not contain harmonics of n = (0, 1), no such restrictions apply 
to the disturbances to be input to Hotine's integral formula. It should be noted here, that 
the absence of first-degree harmonic from T, which implies coincidence of the center of 
mass of the Earth with that of the reference ellipsoid, was implicitly assumed when the 
centrifugal terms of the true geopotential and the normal one were equated (see equations 
(2.8) and (2.10)). The effect however, on the centrifugal potential, of non-geocentricity 
of the coordinate system used, is only about 1 x 10' 3 kgalm (for present accuracies on 
the determination of the geocenter), hence the geocentricity condition there, only mildly 
needs to be employed. 


2.3 Consideration of Systematic Effects 

The integral formula (2.64) was developed based on the assumptions: 

(a) The disturbing potential T is harmonic outside the sphere (O, R) (see equation 
(2.52)). 

(b) The boundary values 8g represent radial derivatives of T (with opposite sign; see 
equation (2.35)). 

(c) Boundary values are given on the surface of the sphere (O, R). 

In the real world, the presence of the atmosphere violates assumption (a), while the 
available boundary data do not comply with either (b) or (c), as was discussed in section 

2.1. The simplicity of the previous formulation though, and the magnitude and spectral 
content of the discrepancies between the real world case and the one assumed in section 

2.2, suggest that it is preferable to retain the formulation developed, and modify the 
observable, in order to ensure the best possible compatibility with the analytical models, 
at the cost of degrading somewhat the integrity of the original measurements. 

The primary information available from surface gravimetry combined with GPS 
positioning, is the magnitude of the gravity acceleration at the surface point P, and the 
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geocentric Cartesian coordinates of P. Normal gravity at the same point may be 
calculated by transforming the Cartesian coordinates to ellipsoidal ones (Heiskanen and 
Moritz, 1967, section 1-19), and then making use of closed formulas for the magnitude 
of normal gravity, as described in (ibid, section 6-2). Alternatively, Cartesian 
coordinates may be transformed to geodetic ones, and normal gravity may be calculated 
from a truncated series (ibid, section 2-10). Pavlis (1988, section 3.1.1) has shown that 
such series should include terms at least up to 0(h^/a^), to avoid introduction of 
undesirable systematic errors. That is: 

|?p| = y Q( [l -2(l +f + m-2fsin 2 (p Qo )^+3{ 1 r) 2 ] - (2.76) 

where <poo is the geodetic latitude and hp the ellipsoidal height of point P (see also Figure 
2), and the rest of the notation is defined in (ibid, section 2.3.1). 

According to the above the gravity disturbance 5gp, defined in equation (2.33), is 
evaluated. The systematic corrections to this quantity are described next. 

1. Atmospheric Correction : 5g A 

Provided the mass of the reference ellipsoid used to define the normal potential 
includes the total mass of the atmosphere (as in the cases of GRS '67 or GRS ’80), 
Moritz (1980, p. 424) has shown that the atmospheric correction on measured gravity is 
given by: 

8 gA = G^ (2.77) 

r p 

where M(r P ) is the mass of the atmosphere above a sphere passing through the 
observation point P (see also (Pavlis, 1988, Figure 6)). To remove the atmospheric 
effect from 8g P , Sg A needs to be added to the gravity disturbance as given by (2.33). 
Obviously, since 8g A has to do with measured gravity, which enters in the same way in 
the definition of both disturbances and anomalies, the correction is identical for both 
quantities (Moritz, 1974). 

For computer implementation it is convenient to evaluate dg A from: 

Sg A = 0.8658 - 9.727 x 1Q- 5 H P + 3.482 x 10 9 hJ (mgal) 


(2.78) 
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where Hp is the orthometric height (or equally well the ellipsoidal one) of the gravity 
station in meters. This polynomial was developed by Wichiencharoen (1982), by fitting 
a quadratic function to the tabulated values of (IAG, 1971, p. 72). The indirect effect 
(shifting of equipotential surfaces due to the condensation of the atmospheric mass on the 
reference ellipsoid), is only about - 0.7 cm at sea level, where the correction becomes 
maximum (Moritz, 1980, p. 425), and thus it can be safely neglected. 

2. Ellipsoidal Corrections : eh, e p 

From equation (2.29) it can be seen that: 

_/3W\ IdW) 

£p \ 3h jp \dH/p (2.79) 

i.e., ep arises due to the difference in the directional derivatives of W along the straight 
ellipsoidal normal (which almost coincides with the normal plumb line) and the actual 
plumb line. On the other hand, from (2.26) it is seen that eh represents the difference (up 
to 0(e 2 )) between the directional derivatives of T along the radial line and the straight 
ellipsoidal normal. Pavlis (1988, sections 2.3.3, 2.3.6) has shown that the corrections 
Eh and ep are almost identical both in terms of magnitude and frequency content. The 
corrections are of the order of 10 pgals and produce long- wavelength signatures on the 
disturbing potential that may reach 20 kgalcm (ibid. Figure 42). The combined 
correction (eh + ep) can be computed efficiently in terms of either point values or area- 
mean values, from an existing geopotential coefficient set such as the OSU89B (Rapp 
and Pavlis, 1990), using the formulation developed by Pavlis (1988, sections 2.3.3, 
2.3.4, 2.3.6). 

With the above correction terms defined, making use of equations (2.30) and 
(2.52), one can write: 


LJk e, a.) = gm. £ (n + iWr £ c^y Je, x) 

Ip 11=0 ' “ 


(2.80) 


where the reduced observable L a (rp, 0, A.) is given by: 
L a (rp, 0, X) = Igpl - lyH + (figA - Eh - £p) • 


(2.81) 
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Equation (2.80) represents a mathematical model, linear with respect to the 
coefficients C nm . A corresponding form of (2.80), written in terms of area-mean values 
of L a (rp, 9, X), may be used to set up a linear system of observation equations, in order 
to estimate a truncated set of coefficients Cnm, from a (preferably global) set of 
"observed" values of La(rp, 0, A,). This may be done in essentially the same manner as it 
was done, using gravity anomalies, in (Pavlis, 1988). 

However, for the current purpose, the reduced observable of (2.81) is clearly 
inadequate, since it refers to the surface point P, while the intention is to make use of 
(2.64), which requires values referring to the sphere (O, R). The continuation of the 
values referring to the physical surface of the Earth, to a surface that is analytically 
manageable, poses the most difficult problem of all other reductions, both from a 
theoretical as well as a numerical standpoint. The treatment described next consists of 
two steps; first the surface values are analytically continued to the surface of an ellipsoid 
and second, equation (2.67) is modified to account for the differences between the 
ellipsoidal and the spherical surfaces. The problem of analytical continuation is also 
encountered in the implementation of Stokes' integral, and due to its importance has 
received extensive studying by a number of investigators (e.g. Moritz (1966; 1974; 
1980); Wang (1988)). Their established notation is adopted in the following discussion. 

3. Analytical Continuation : gj 

For notational convenience the quantity L a (rp, 0, A) of (2.81) will be simply 
denoted 8g here (not to be confused with 8g as used in other sections). The purpose of 
analytical continuation, is to determine a corresponding quantity 8g* on the ellipsoid, 
such that, 8g is related to 8g* through Poisson's upward continuation integral (Heiskanen 
and Moritz, 1967, p. 35). The quantity 8g* possesses no physical meaning, as the 
surface value 8g does, its only purpose being to enable use of convenient analytical 
formulas. Relating the two quantities through the upward continuation integral, ensures 
that the set of 8g* values on the ellipsoid, reproduces the set of 8g values on the surface, 
and consequently the external potential determined from 8g* is identical to the one that 
would have been determined from 8g. Strictly speaking, the use of Poisson's integral is 
valid only if the reference surface of 8g* is always below the physical surface of the 
Earth, a condition which is not satisfied (in general) by the surface of the mean-Earth 
ellipsoid. This problem led Bjerhammar (1962) to introduce the concept of a spherical 
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reference surface for 8g*, completely embedded in the Earth. The validity of Poisson's 
integral is also ensured if one introduces an embedded ellipsoid, confocal with the mean- 
Earth one, as a reference surface for 8g*. This way analytical continuation has to be 
performed over shorter distances, which is advantageous from the accuracy and 
convergence standpoints, while the embedded ellipsoid concept does not add significant 
complexity to the formulation; all that it requires is to reckon ellipsoidal heights from the 
embedded, instead of the mean-Earth ellipsoid in the following formulas. 


Let now Rb denote the radius of an embedded sphere. Application of Poisson's 
integral to the harmonic function r5g yields: 



(2.82) 


where the notation could be read from Figure 3 (substituting Rb for R used there). Since 
5 g * is the unknown quantity, the integral equation (2.82) needs to be inverted. 
However, there exists no integral formula that inverts (2.82) (Heiskanen and Moritz, 
1967, section 8-10), hence the solution of (2.82) can only be obtained numerically, with 
successive approximations. To this end, it can be shown easily (ibid, p. 318) that (2.82) 
may be written in the form: 


K 2j. * t^l- 1 2 ) 

5gp-t8g P =-^ 


8g - s g ; 

D 3 


do 


where : 


t = Rfe. 
rp 


D = f 

rp 


(2.83) 

(2.84) 


The location of the '*' quantities on the embedded sphere is determined by projecting the 
surface points to this sphere along the radial line, i.e. the surface point and its projection 
have the same geocentric latitude. Approximating 


Tp ~ Rb + hp , 


(2.85) 
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expanding t 2 and D 3 in powers of hp/Rb, and retaining only terms linear in h P /Rb, one 
can easily show that (2.83) takes the form: 


5 g ;> 


5gp - hp - =2- 8g P + f f — , § — da 


R b 


io 


where 


( 2 . 86 ) 


to 


= 2R b sin y 


(2.87) 


Equation (2.86) lends itself to an iterative scheme for the computation of 8gp which 
is initialized by setting 8gp<°) equal to 8g P , as described for the case of gravity anomalies 
by Heiskanen and Moritz (1967, p. 318). Moritz (1966, p. 60) has shown that under the 
assumptions made above to derive (2.86), and to the first order of hp/Rb, this equation 
coincides with the 'gradient solution' to the analytical downward continuation problem, 

i.e. the terms inside the brackets in (2.86) represent (^8 /dr). The numerical 
implementation of (2.86) on the other hand is all but trivial, since the vertical gradient 
required has a veiy localized behavior and thus its accurate estimation requires very dense 
observation coverage around the computation point P. Usually, such coverage is not 
available. Employing assumptions regarding the correlation of the observable with 
elevation, and making use of available detailed elevation information, an approximate 
solution to (2.86) may then be evaluated, as was done by Wang (1988) for the case of 
gravity anomalies. The problem of analytical downward continuation is a topic on its 
own and further discussion will not be given here. Equation (2.86) will be written 
schematically: 

8gp = 8gp + gi (2.88) 

where : 


gi = - hpj 




5g*- Sgp 

4 3 


da 


(2.89) 
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As far as the use of an embedded ellipsoid is concerned, Rb in (2.89) may be 
substituted by the Gaussian mean radius (Rapp, 1984, p. 43) at the computation point P, 
for that ellipsoid, since global integration in (2.89) is usually truncated only over a small 
cap centered at P. 

4. Ellipsoidal to Spherical Integration : ^ 


The values 8g* obtained from the previous step refer to the surface of the 
embedded ellipsoid and thus are still inadequate to be used as input to Hotine's integral 
(2.67), which requires values on a spherical boundary surface. This problem, for the 
corresponding case of Stokes' integral, is usually treated with ellipsoidal correction terms 
as those derived by Moritz (1980, p. 314). However, as Hotine (1969) has pointed out, 
the problem may be treated in a conceptually simpler manner; his formulation is 
presented next 


First, it is observed that the result of downward continuation, 5g*, represents 
radial d erivative of the analytical continuation of the disturbing potential (see also 
equation (2.81)). Ellipsoidal normal derivatives will be needed next, which formally may 
be obtained by adding back (eh)p to 5gp, or, without any loss of accuracy by omitting 
altogether (eh)p from (2.81). Hence, 5gp from (2.88) is hereon understood to represent: 



With respect to the ellipsoidal coordinate system (a = cot' lu /E, 5, A.) (see also (Heiskanen 
and Moritz, 1967, sections 1-19, 1-20)), the ellipsoidal harmonic expansion of the 
disturbing potential will be written as: 




where: 

i=VT 


(2.91) 


cota b = 4- 


(2.92) 
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and eb is the second eccentricity of the embedded ellipsoid, while Qnm( z ) are associated 
Legendre functions of the second kind with argument z (ibid, p. 43). ¥^(5, X) are as in 
(2.53), but now evaluated in terms of reduced co-latitude 5. Finally are real , 
ellipsoidal harmonic coefficients of the disturbing potential, referring to scaling 
parameters GM and Rb. Making use of the relation (Hotine, 1969, p. 190): 

JL , . tan g JL 

3h N da (2.93) 


where N is the radius of curvature of the prime vertical (Rapp, 1984, p. 32), and the 
recurrence relation for the derivative of Qn m ( z )» one can derive easily: 


N aT _ y y Tnmfs, X) 
a h m-n Qnln^icotCtb) 


(n + ljQjjn/icota) + itanot(n - m + l)Q,+i (icota) 


(2.94) 


where the surface ellipsoidal harmonic Tnm(S, X) is: 

T nm (8,x) = ^C e nm Yj8,X) . 


(2.95) 


Equation (2.94) holds true for any arbitrary point on, or outside the embedded ellipsoid. 
On the surface of this ellipsoid, one has N = Nb and a = <Xb, hence due also to (2.90), 
(2.94) becomes: 


oo n 


Nb8g = £ X ( n + l)Tnm(8, X) 

[itanab(n - m + l)Qi+i m (icototb) 


n=o m=-n 
oo n 

1 1 

n=o m=-n l 


Qi|m|(icotab) 


Tnm(8, X) 


(2.96) 


Using the series expression for Qnm( z )» it can be shown that the bracketed term in the 
second summation of (2.96) is equal to: 


e^ 2 (n - m + lXn + m + 1) 
2n + 3 



(2.97) 
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The last relation is a corrected version of the misprinted relation (29.66) of Hotine (1969, 
p. 321). Hence, omitting terms of the fourth and higher order of the second eccentricity, 
equation (2.96) may be written as: 


oo n 


[N b 8g*](8', X') = X £(n + X') 

n=o m=-n 

+ c' b 2 £ £ ( n - m+ 'X" + ™ + i> Tjs-. xj 

11=0 m=-n in + 5 


(2.98) 


Multiplying both sides of the last equation by H(y) as given in (2.65) kul with \|/ now 
evaluated by: 


cost)/ = cosScosS' + sin5sin5'cos(A. - X') , 


(2.99) 


and integrating over the unit sphere (full solid angle), due to (2.66) one obtains: 


l(6, X) = ±jj H(vfN b Sg-] (S', X'jdo 


4 K I | n=o m=-n 


£ £ (n - m ^; 3 +m+1) h( V )tJ 8'. r)dc 


Denoting: 


a(8', X') = e' b 2 £ £ ( - m V X , n 3 m + 1 > T J8', X’) 

n=om=-n zn + ^ 


and collecting all the previous correction terms together, one finally obtains: 


( 2 . 100 ) 


( 2 . 101 ) 


I Hw 


N b (Sg + SgA - ep + gi) - a 


(s', X')dc 


( 2 . 102 ) 
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Equation (2.102) is the desired integral formula which relates the gravity 
disturbance 5g (as defined in equation (2.33)), to the analytical continuation of the 
disturbing potential on the surface of the embedded ellipsoid, accounting for the ellipticity 
of the boundary surface to the second order of the second eccentricity. It can be seen 
now that the use of spherical formulas (such as (2.67)), but evaluated in terms of reduced 
instead of geocentric colatitudes, introduces two errors; one due to the difference 
between Nb and R, and another due to the omission of e,. Both errors are of the second 
order of the (second) eccentricity (Hotine, 1969, p. 321). Also, from (2.101) and (2.95) 
it is obvious that the computation of e, requires some a-priori knowledge of the 
ellipsoidal spectrum of the geopotential. To avoid the use of the ellipsoidal harmonics, 
and subsequently the need for transformation of geodetic latitudes (according to which 
measurements are usually registered), to reduced ones, Moritz (1980) used series 
expansions and formulated the ellipsoidal corrections, for various gravimetric quantities, 
in terms of geodetic latitudes and spherical harmonics (ibid, pp. 314-328). Such 
procedure, as Moritz himself admitted, results in more complicated expressions, than 
those for ellipsoidal harmonics (ibid, p. 316). At present, neither one of Moritz's 
motives appears to be a prohibitive factor for the formulation adopted here. The 
derivations of Jekeli (1988) and the computer algorithms of Gleason (1988), provide 
very efficient means of converting ellipsoidal to spherical spectra and vice-versa. 
Furthermore, in view of the computational facilities available nowadays, conversion of 
geodetic to reduced latitudes, even for large amounts of observation locations, can hardly 
be considered a prohibitive computational task. 

In practical applications the integration in (2.102) is usually extended only over a 
small cap, centered at the computation point, while one accounts for the effect of the 
remote zone through the use of a global geopotential model. Such procedure requires 
modifications in equation (2.102), which affect the evaluation of e, as well, as it will be 
discussed in section 3. 1 . However, it is instructive here to evaluate the effect of ^ on T, 
for the case of global integration. 

Denote: 

tf(s, II H|v)s(5'. v)d<j 


(2.103) 
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Substituting in equation (2.103) £i( 5', X') from (2.101) and H(\y) from (2.66), and 
making use of the decomposition formula of the Legendre polynomials (Heiskanen and 
Moritz, section 1-15), one obtains: 


ef(S,X) = < 2 £ £ 

n=0 m=-n 


(n-m+1) (n+m+1) 
(2n+3) (n+1) 


Tnm(5,X) . 


(2.104) 


Given now a set of spherical harmonic coefficients of the disturbing potential, Cnm. 
referring to scaling parameter R = R M , one sets: 


f* 

' nm 


= GM 
Rb 



(2.105) 


Using the transformation formula (1.15) of Gleason (1988, p.116), one can evaluate 
from /nm, the ellipsoidal spectrum while the unitless ellipsoidal coeffecients C£m are 
finally given by: 

c e - Rb 7 e 

Wn GM Jnm ' (2.106) 


The term ef(5,X) may now be evaluated from equations (2.95) and (2.104), using 
efficient harmonic synthesis techniques such as those developed by Colombo (1981a). 

Such evaluation was performed here using the OSU89B shperical harmonic 
coefficient set (Rapp and Pavlis, 1990), complete to degree and order 180; ef(8,X) was 
evaluated in terms of 1° x 1° area-mean values, while for simplicity Rb was set equal to 
R M = 6378137 m, the scaling parameter to which the OSU89B model coefficients refer 
(ibid, p. 21,896). The ellipsoidal gravity field to which Cnm refer was defined through 
the four constants given in (ibid, p. 21,896). The values of ef(8,X) computed, range 
between -19 kgalcm and +17 kgalcm, with a root mean square (rms) value of 6.4 
kgal cm. Their geographical variation is illustrated in Figure 5. 
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Figure 5. Disturbing Potential Correction ef. Area-Mean Values (1° x 1°) Computed Using OSU89B to Nmax = 180. 
Contour Interval is 2 kgal cm. 
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In a similar fashion one can derive the corresponding correction term for the case 
of integration of gravity anomalies (Stokes' integral); such derivation yields (see also 
(Hotine, 1969, p. 322)): 


(^)a 8 (5, X) = e' b 2 X 2 

n=0 m=-n 


(n-m+1) (n+m+1) 
(2n+3) (n-1) 




oo 


I 


X ( e2 sin 2 S + T nm (5,A.) 

n^n n - 1 \ Y I 


(2.107) 


where the prime indicates absence of the first-degree term from the infinite sum. 
Numerical values computed for the correction term of equation (2.107) (in the same 
manner as for the term of equation (2.104)) range between - 36 kgal cm and + 43 kgal cm, 
with an rms value of 14.9 kgal cm. Their geographical variation is illustrated in Figure 6. 
These values are in good agreement with the corresponding ones computed by Rapp 
(1981b), using the procedures of Moritz (1980). Both correction terms (2.104) and 
(2.107) are of long- wavelength nature, with more than 99% of their power concentrated 
below harmonic degree 36, so that one can account for their effect accurately, using an 
existing global geopotential model. 


Finally, to evaluate the disturbing potential at the surface point P, one needs to 
upward continue the value computed from (2.102) (which refers to the footpoint P* on 
the embedded ellipsoid), to the topographic surface level. A Taylor series expansion in 
terms of the elevation yields: 


T(P) = T(P*) + 




(2.108) 


hence, due to (2.90) and (2.89) one finally has: 
T(P) * T(P*) - (8gp - ^ gi)h P . 


(2.109) 
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Figure 6. Disturbing Potential Correction (e")Ag. Area-Mean Values (1° x 1°) Computed Using OSU89B to Nmax = 180. 
Contour Interval is 2 kgal cm 
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2.4 Gravitational Acceleration Information from a Satellite-to-Satellite 

Tracking Configuration 

The other source of gravitational information considered here for the estimation of 
the geopotential differences, comes from the Multiple-High-Single-Low Satellite-to- 
Satellite Tracking (SST) configuration which is established when a low orbiter carrying a 
GPS receiver is simultaneously tracking more than one satellite of the high-altitude GPS 
constellation. 

The idea of using SST data for geopotential modeling was originally proposed by 
Wolff (1969). The fact that an SST low- low mission alone, with the satellites in polar 
orbit, would be capable of providing a truly global, uniformly accurate and high 
resolution geopotential model, caused widespread attention to be given to the proposal. 
A number of investigations have been performed, aiming to assess the capabilities of 
various SST system scenarios. Some of these studies aimed to assess the quality of 
mean gravity anomalies and/or geoidal undulations that can be predicted on the surface of 
the Earth from the SST data at altitude (local solutions) (e.g. Hajela( 1974; 1978; 1981), 
Rummel et al. (1976), Rummel (1980), Rapp and Hajela (1979) and Douglas et al. 
(1980)). Other studies (e.g. Colombo (1981b), Kaula (1983)) focused on the 
development of efficient analytical techniques, that can be used to process the large global 
set of observations which an SST mission will acquire during its lifetime, for the 
determination of harmonic coefficients of the geopotential. However, apart from early 
experiments of SST between ATS-6/Apollo-Soyuz and ATS-6/GEOS-3 (Kahn et al., 
1982), no dedicated SST system has yet been put into orbit. This is partly due to the 
costly requirement for drag-free satellites at the low altitudes (e.g. 160 km), considered 
for SST missions of the type of the Geopotential Research Mission (GRM) (Keating et 
al., 1986). 

In the absence of a dedicated SST mission, the possibility of using the GPS 
constellation as the high-altitude component, which continuously tracks a low-altitude 
satellite equipped with a GPS receiver, offers a viable alternative. Investigations in this 
direction have been already initiated through the U.S. Air Force "Shuttle-GPS Tracking 
for Anomalous Gravitation Estimation" (STAGE) mission (Jekeli and Upadhyay, 1990), 
where the Shuttle spacecraft is used as the low-orbiter, and an Inerdal Measurement Unit 
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(IMU) is used in addition to the GPS receiver on the Shuttle, to measure and isolate non- 
gravitational accelerations. 

The basic mathematical modeling of the observations acquired by an SST 
configuration remains unaltered regardless of the mission in question; in contrast, the 
signals represented in the observable, strongly depend on the specifications of each 
mission (e.g. drag-free orbit, altitude etc.). The modeling of the observable is reviewed 
next, while the contribution of various signals contained in it is discussed afterwards. 

Consider the motion of two satellites, Si (high) and S 0 Cow), as shown in Figure 
7. Three mutually perpendicular unit vectors Ej (j = 1, 2, 3) span an inertial frame of 
reference (in this section, vectors will be denoted with underbars, for notational 
convenience). 



Figure 7. High-Low Satellite-to-Satellite Tracking Configuration. 


The following relations arc introduced: 
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|liol = (lLlio) l/2 = Pio 


Cio = 



^io 


e T io nio=0 . 

Denoting with d/dt time differentiation in the inertial frame, one has: 
iio = ^(-Tio) » 


( 2 . 110 ) 

( 2 . 111 ) 

( 2 . 112 ) 


(2.113) 


for the inertial relative velocity between S 0 and Sj, and : 

^io = Pio — Eio-fiio » (2.114) 

for the projection of the inertial relative velocity along the direction of the "line-of-sight" 
(LOS) between the two spacecrafts. Designated "LOS velocity", v io is a quantity that can 
be observed from the Doppler effect on the signal received by the low orbiter. 
Differentiation of (2. 1 14), with respect to time yields: 

Vio = f io£io+ tio £ io • (2.115) 

The quantity: 


aio ** Xio £io (2.116) 

is the projection of the inertial relative acceleration along the LOS direction. On the other 
hand,Vio is the rate of change of the LOS velocity. The second term on the right-hand 
side of equation (2.1 15) (centrifugal acceleration) arises due to the fact that e io represents 
a direction which is rotating with respect to inertial space. 
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Assume now that a set of initial state-vectors for both satellites and force model 
parameters for the forces acting on them are approximately known. Based on such 
information one can compute an approximate value for Vi c as: 

v? 0 = (i1o) T ef 0 + (ifo) T fi io • ( 2 . 1 17) 


However, since neither the initial state-vectors, nor the force model parameters are 
known perfectly, one has: 


5vj 0 = Vj 0 - vf 0 , or 
Svio = 6(.Tio£io) 8(iio£ ic) 


where: 


S(f io£io) — I ioSio " (l io) £* 


,c 

io 


8(f ioSio) — EioSio " (iio) £ 


,c 

io 


(2.118) 

(2.119) 


(2.120) 


One may now split the total inertial relative acceleration jr io into two parts: 

~g . ..no 

lio ~ Xio + lio > (2.121) 

where i}° is induced by gravitation, and r^° is due to all other forces Jim gravitation, 
acting on the spacecrafts Sj and S 0 . Accordingly, equation (2.1 19) takes the form: 

Svio = 5[(f°o) T £io] + 8[(f^) T £io] + 5 [i ioSio] • (2.122) 


The first term on the right-hand side of (2.122) will be the focus of the following 
discussion. As mentioned before, this term arises due to the errors of the approximately 
known gravitational model and satellite positions. Assuming that these errors are small 
enough to justify linearization, one has: 


S(fIo£io) * (§i 


10 




^(f iofiio) 


3r, 


'io 


5io + 


afel 


o-£io; 


3ii 


5n 


C 


c 


(2.123) 
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where the superscript "G" has been omitted for notational simplicity, and in the following 
it will be understood that refers to the component of the relative inertial acceleration 
induced by gravitation alone. In addition "c" indicates that the partial derivatives are 
evaluated at the approximately known positions. 


On the other hand, by differentiating numerically the observed Doppler shifts, one 
obtains an "observed" value for v io , v °o s , where: 


v?o = Vj 0 + n io (2.124) 

and njo represents observational noise in the primary observable, unmodeled (residual) 
atmospheric effects, antenna multipath, as well as errors introduced by the numerical 
differentiation. Introducing the notation: 


Avio = v» s - Vio ( a ) 

5a^ = (i^) T £io-[(iOTfif 0 (b) 

5R io = 5(tioeio) (c) . 


(2.125) 


one has from (2.122), (2.123) and (2.124): 


Av 


10 


= (5iio) 1 


£io + 


^(x io-fijo) 


die 


5lo + 


3(r 


Uio) 


dli 


5ri+5a^ + 6R i0 


+ e 


10 


(2.126) 


where £i 0 contains now apart of the noise n; 0 , the effects of second and higher order 
terms omitted in equation (2. 123). For the first partial derivative in (2.126), one has: 

3(.Tio-£io) _ _T 3(Xio) , i;T afeiol 

Bio _ i ° dlo Iio d_r o • (2.127) 


It can be shown easily (see also (Rummel, 1980)) that: 
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where I3 is the three-dimensional unit matrix. Rummel (ibid) also shows that: 

dlo dr o (2.129) 

where V is the (true) gravitational potential of the Earth, and M the three-dimensional 
gravitational tensor, i.e., 


M(r c ) = 


8 2 V 

d 2 V 

d 2 V 

dx 2 

dxdy 

dxdz 

8 2 V 

d 2 V 

d 2 V 

dydx 

dy 2 

dydz 

a 2 v 

d 2 V 

d 2 V 

dzdx 

dzdy 

dz 2 


(2.130) 


With corresponding derivations for the partial derivative with respect to rj, equation 
(2.126) after regrouping terms becomes: 

Avi 0 = feojsfo + (Ufo) T M(lg) + ^-[i:‘Io - aiotefofljsio 

Cef 0 ) T M(i?) + -j- [ilo - aioCefofljsii 
Pio / 

+ 8a? 0 ° + 5R i0 + eio • (2.131) 

From equation (2.131) it can be seen that the (pseudo) observable Av io is affected 
by (and thus contains information concerning) the following: 

(a) The residual (with respect to the reference geopotential model used) relative 
gravitational acceleration along the LOS between the two spacecrafts (first term). 


(b) The difference between the actual and the reference orbits of the two spacecrafts 
(second and third terms). 
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(c) The residual (with respect to the non-gravitational force models used) relative non- 
gravitational acceleration along the LOS, 8a^. 

(d) The residual relative centrifugal acceleration along the LOS, 8Rj 0 . 

(e) The random noise on the primary observable, residual atmospheric effects, non- 
linear term effects etc., all these effects collectively represented in the term e^. 

For the current application, of interest is only (part of) the first term in equation 

(2.131) , while the rest of the terms represent undesired systematic and random "errors". 
In contrast, if for example the intention was to estimate improved orbits for the two 
spacecrafts, based on the (pseudo) observations Av io , then the second and third terms in 

(2.131) would have been of critical importance. Since: 

5i i0 = 5i 0 - 8ii , (2.132) 

and: 

Sfiio = -L-[h- (fifo) (efo)1 (5io - 5ii) , 

Pio (2.133) 

one can re-write equation (2.131) as: 

Avjo = (5r'o)Vo * (Sfiffifo 

+ (ef 0 ) T [M(r8)8r 0 - M(r?)5ii] + r[ 0 8eio 

+ 8a?° + 8R i0 + £io , (2.134) 

where the orbit error contributions were separated in the parts referring to the absolute 
and relative position errors of the spacecrafts respectively. The magnitude of each term 
on the right-hand side of (2.134) depends on the particular SST configuration in question 
(e.g. satellite altitudes, instrumentation etc.), the maximum degree (as well as the 
accuracy) of the reference geopotential model used to evaluate Avi 0 , and the accuracy of 
the satellites’ ephemerides. The following considerations pertain to the case where S 0 
corresponds to GP-B (see Table 1), while Si corresponds to one of the satellites of the 
GPS constellation, so that: 

r 0 =|lol * 6971 km 


r i = llil “ 26560 km 
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To estimate the magnitude of each term in (2.134) one needs to assume some knowledge 
of the terrestrial gravitational field, at least in a global average sense, as described by a 
global covariance model (Moritz, 1980, p. 181). The model used here will be defined by 
the following anomaly degree variances (see also section 3.4): 


c„(mgal 2 ) = ( 






n+2 


2 < n < Nmax 


Nmax < n < 360 


360 < n < oo 

(2.135) 


where: 


GM = 3986004.36 x 10» 

m 3 /s 2 

a = 6378137. 

m 

R = 6371000. 

m 


and the C nm , eCnm are the fully-normalized unitless harmonic coefficients and their errors 
as given by the OSU89B geopotential model (Rapp and Pavlis, 1990) (even zonal 
harmonic coefficients are remainders after removing the ellipsoidal reference field). The 
c„ values defined above refer to the surface of the sphere of radius R, while Nmax 
represents the maximum degree of the geopotential model used to evaluate Avj 0 . 

Considering now the first term on the right-hand side of equation (2.134) one has: 

|(5r 0 ) T £f 0 | <|Sf 0 | = |grad[T m (r8)]| = (2.136) 

where T™ is the disturbing potential with respect to the reference model used (i.e., T™ 
represents the commission error of the model up to degree n = Nmax, and the omission 
error from n = Nmax + 1 to infinity), and 5g m is the gravity disturbance implied by the 
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aforementioned model (i.e. 5g m = - ^ /dr )• The difference in direction between the 
actual gravitational acceleration and the one implied by the model has been disregarded in 

(2.136). In a spherical approximation, which is sufficient for the magnitude estimates 
sought here, one has (see also (Jekeli, 1979)): 


rms(o) = (var[ (SiofifJ } ‘ 



(2.137) 


and similarily, for the high-altitude satellite: 


rms(i) ■ (var[ (6f iffifj f* m 



(2.138) 


The cross-rms of the residual accelerations of the high and low satellite, assuming radial 
arrangement of the spacecrafts, as the \|/ = 0 notation indicates, is given by (Jekeli, ibid): 




(2.139) 


This cross-rms is maximized for radial arrangement of the satellites, enabling a worst 
case study of the effect of its omission. In addition, the rms acceleration difference 
between the low and high satellite (assuming radial arrangement), is given by: 

rms(o - i) = (var[ (Si'io * [rms^o) - 2rms 2 (i, o) + rms 2 (i)]^ . (2.140) 


These quantities have been evaluated using the degree variances defined in (2.135), and 
approximating infinity by n = 36000 in (2.137) through (2.139). The results for various 
degrees of truncation (the variable Nmax appearing in equation (2.135)), are given in 
Table 2. 
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Table 2. Root Mean Square (rms) Residual Acceleration Magnitude With Respect to 
OSU89B Model Complete to Degree Nmax (All rms values are in mgals). 


Nmax 

rms(i) 

rms(o) 

rms(i, o) 

rms(o - i) 

2 

0.94106 x 10- 2 

10.748 

0.27752 

10.741 

4 

0.32439 x 10-3 

6.391 

0.37761 x 10 1 

6.391 

6 

0.16776 x 10- 4 

4.228 

0.68755 x lO- 2 

4.228 

8 

0.55717 x 10-5 

2.784 

0.11287 x 10-2 

2.784 

10 

0.55323 x 10- 5 

1.865 

0.25036 x lO’ 3 

1.865 

20 

0.55323 x 10-5 

0.583 

0.18081 x lO 3 

0.583 


From Table 2 it can be seen that a state-of-the-art reference gravitational model 
(developed in the absence of any of the missions discussed in section 1.3) and truncated 
to a degree as low as eight is enough to justify the assumption that the residual 
gravitational acceleration at GPS altitude is zero, thus introducing an error no larger than 
about 10-3 mgal (see also (Jekeli and Upadhyay, 1990)). It should be noted that Jekeli 
and Upadhyay (ibid), consider the reference model up to Nmax errorless, thus showing a 
monotonic decrease for rms(i) as Nmax increases; here rms(i) stabilizes at about 
0.6 x 10' 5 mgal after Nmax = 8, due to the commission error of the lower degree 
harmonics. 

The effect on Avj 0 of the orbit errors 8r 0 and Sri of the low and high satellites is 
considered next. For this purpose, it is assumed that a global network of tracking sites 
on the Earth, simultaneously co-observes the GPS satellites being tracked by the low 
orbiter, as it will be the case for all the missions discussed in section 1 .3 (Pavlis, E. et 
al., 1990). In such case the orbits of the GPS satellites and the low orbiter can be 
estimated geometrically from differential GPS observations. Yunck and Wu (1986) 
carried out simulation studies to assess the accuracy of such orbit determination for 
TOPEX, and concluded that decimeter accuracy orbits are attainable with as few as 10 
globally distributed tracking stations. The orbital accuracy in such non-dynamic 
solutions is limited by the GPS data noise, the ground station position errors, 
tropospheric and higher order ionospheric effects and antenna multipath, error sources 
that have no significant dependence on altitude (ibid, 1986), so that 10 cm is also a 
realistic (and probably too conservative) estimate for the orbit error of GP-B. 



Approximating the gravitational tensor M(i) by: 
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Mi (l) 


_£M 

Id 3 


10 0 ' 
0 1 0 
.0 0 - 2 . 


(2.141) 


which is sufficient for order of magnitude considerations (see also (Jekeli and Upadhyay, 
1990)), and denoting by e r the radial component (which is the most crucial for the current 
application) of the acceleration error induced by radial orbit error 5r, one has (see 
equation 2.134): 

Low satellite (S 0 ) : (mgal) * 2.4 x lO^Sr^m) 

High satellite (Si) : e[ (mgal) » 4.3 x 10' 3 8rj(m) 


while the misregistration of the LOS direction, 8£ic introduces an error £ 5 ^ > which 
does not exceed: 

e§ eiD (mgal) = 4.2 x 10' 2 V2 8ro(m) 

As it can be seen from the above estimates (which pertain to the GPS/GP-B SST 
configuration), the orbit error of the GPS satellite and the misregistration of the LOS 
direction introduce negligible acceleration errors. The orbit error of the low satellite 
however, introduces an acceleration error that may reach 24 (igals in amplitude (for 8r 0 = 
10 cm) which can have a significant effect on geopotential estimation at the few 
centimeter accuracy level, if its power is concentrated at low frequencies. Colombo 
(1990) has shown that systematic errors in the satellite positions significantly affect the 
recovered geopotential spectrum only up to about degree 10. In addition, comparing the 
results of the analytical approach in the absence of such errors ("best case"), with those 
obtained from a complete simulation where satellite orbits were estimated dynamically in 
a simultaneous solution along with the geopotential spectrum (Pavlis, E. et al., 1990), 
Colombo (1990) verified that systematic orbit errors can be effectively decoupled from 
the gravitational signal in such global solutions, as long as their mathematical 
representation is accounted for in the adjustment. Accordingly, one way of reducing the 
effect of orbital errors on Avj 0 , is to pre-process smoothed values of the global set of 
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SST data in a dynamic mode where the orbital parameters of all satellites and the low- 
degree part of the geopotential spectrum (e.g. Nmax = 10) are simultaneously estimated 
from the original GPS measurements. Avj 0 can then be referred to the adjusted orbits 
and long-wavelength geopotential spectrum obtained from such global dynamic solution. 
In such case, the smoothing of the original measurements is critical, since one wants to 
minimize leakage of the higher-frequency content of the measurements to the lower- 
frequency part of the estimated geopotential spectrum. However, the global character 
(polar GP-B orbit) and the uniform coverage of the data works favorably in the 
minimization of leakage effects. Obviously, as the maximum degree of the geopotential 
model obtained from the global adjustment increases, the contribution of the residual Avj 0 
(with respect to this model), to the estimation of the disturbing potential T on the Earth's 
surface, decreases. In the limit, if the model extends to the degree corresponding to the 
resolution of the data at altitude (Nmax « 55 for GP-B), then only the ground 
measurements, in the caps surrounding the computation points provide additional 
information for the estimation of T. An alternative way of reducing the orbit error of the 
low satellite is the use of laser ranging combined with the GPS tracking in a geometric 
solution for its orbit determination (Everitt et al., 1989). In view of the magnitudes of 
8Rio and Ei 0 , which will be discussed next, the effect of orbit errors in the modeling of 
Avj 0 does not appear to be a limiting factor. Accordingly, equation (2.134) becomes: 

Av io « (5r o ) T £f 0 + 8R i0 + £io , (2. 142) 

where, in addition, it was assumed that the residual non-gravitational acceleration along 
the line of sight, Sa^ 5 , is negligible. Such an assumption implies that non-gravitational 
accelerations have been modeled perfectly as a result of the global adjustment of the SST 
data. 


The centrifugal term 8Ri 0 is considered next. It can be easily shown (Colombo, 
1981b), that: 

Rio = ilo §io = -Mi T io n iQ ) 2 = -Mhopcos^ 

Pio Pio 


(2.143) 
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where nio is the unit vector perpendicular to the LOS direction as defined in equation 
(2. 1 1 2), and y is the angle between i io and a io- The geometry of the vectors defined on 

the plane generated by ii 0 and i i 0 is illustrated in Figure 8. 


£io 



Figure 8. Geometry of the Vectors Defined on the Plane of ii 0 and j ;<> 

Subtracting the reference value from the true value Ri 0 one has, due to 
(2.143), for 5R i0 : 

5R i0 = - 1 - |r iofcos 2 ? - -L |rtoP cos V 

Pio Pio (2.144) 

where f is the angle between r f 0 and To obtain a magnitude estimate of 5Ri 0 , it is 
justifiable to set: 

Pio “Pio ; cosy “cosy 0 


so that: 


SR io « 4" { I i id 2 - 1 i foPjcosV = 4- 5{ I i io| 2 )cosY 


'IO 


Pio 


(2.145) 


(2.146) 


8Ri 0 = -2-|rfol|6iio|cos(5cosY > 

Pio 

where the total differential 8(i T i 0 i i 0 ) was approximated by the linear term in its Taylor 
series expansion (linearization), and "(3" is the angle between i f 0 and 5i j 0 . 

The last equation indicates that 8Ri 0 becomes maximum when y = 0°(180°), i.e. 
when the relative velocity of the satellites is perpendicular to the LOS direction. If in 
addition the satellites move in parallel and opposite directions (perpendicular to the LOS), 
then the magnitude of their inertial relative velocity becomes maximum, equal to the sum 
of the magnitudes of their individual inertial velocities. Under such circumstances the 
rate of rotation of the LOS direction is maximized, giving rise to the maximum possible 
centrifugal acceleration Rj 0 , as expected. This "worst case" will be considered next, in 
order to estimate an upper bound for the magnitude of 8Ri 0 ; the "best case" obviously 
corresponds to y = 90°(270°), when the satellites move towards or away from each other, 
so that £io remains fixed with respect to inertial space and Rio = 0. For the "worst case" 
one has from (2.146): 

|SRio| = \ )t c iol|5riol|cosp|cosY 

Pio 

^ IftofSriol 

Pio 

^ -2-(lr c ol +|iil)|8rio| => 

Pio 

max |SRiJ - -2- (| fj + 1 fj) |Sr iol . 

Pio (2.147) 

Making use of the energy conservation law (Jekeli and Rapp, 1980, p. 4), linearized with 
respect to the reference gravitational model to which 8Ri 0 refers, and assuming that the 
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motion of the high-altitude satellite is perfectly determined from that model, so that 
5i i = 0, one has (see also (Rummel, 1980)): 



(2.148) 


where: 


(2.149) 


andT™ carries the same meaning as in equation (2.136) given before. Accordingly, an 
upper bound for the magnitude of 5Ri 0 is given by: 


max|8Rj - -2- 
Pio 



(2.150) 


An estimate of the disturbing potential difference T^, between the low and high satellites, 
may be obtained in a global average sense, considering the radial arrangement of the 
spacecrafts, by: 


1 1/2 

rms(TS) = {var[f , "(r 0 )] - 2 cov[t” (r„), T m ( ri )]^ + var [¥>,)] , (2.151) 


using the anomaly degree variance model of equation (2.135) to evaluate the disturbing 
potential's variances and covariances. Approximating: 



and using the nominal values pertaining to the GPS/GP-B configuration, estimates of 
maxjSRjJ have been computed, for various degrees of truncation, Nmax, of the reference 
model. These are given in Table 3. 
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Table 3. Estimates of the Maximum Value Attainable by 5Rj 0 for Various Degrees of 
Truncation of the Reference Model (Units are mgals). 


Nmax 


5R i0 

■ 

2 



4 



6 



8 



10 



20 

0.023 


1 

1.910 (*) 

18 

0.097 (*) 


(*) Values correspond to the High-Low mission considered by Rummel (1980). 

For comparison purposes, two additional values are listed in Table 3. These refer 
to the High-Low mission considered by Rummel (1980), where the altitudes of the high 
and low satellites were 35500 km and 250 km respectively. These two values have been 
computed using the same anomaly degree variance model as in Rummel's study (ibid, p. 
12). Comparing the above estimate for Nmax = 1 (i.e. when T™ represents the entire 
disturbing potential with respect to the ellipsoidal gravity field) to the corresponding one 
of Rummel's (ibid, p. 5, Table 1), it is seen that the former is about 7 x 10 5 times larger 
than the latter. Colombo (1981b, p. 15) estimated the magnitude of the centrifugal term 
8Rio for a Low-Low mission, and found a value 3 x 10 4 times larger than corresponding 
estimate of Rummel (1980, p. 5, Table 1). 

The reason that both the current estimates, as well as those of Colombo, are so 
much larger than those of Rummel, is identified here to be an unjustifiable substitution of 

the term s[| jr io | 2 ], by the term |5±io| 2 . in Rummel's derivations (ibid, pp. 4-5). Such 
substitution yields an expression for the average magnitude of 8Ri 0 (ibid, p. 5, equation 
12), which is in fact independent of the relative inertial velocity of the two satellites, thus 
contradicting the underlying physical meaning of 8Ri 0 . Also apart from underestimating 
the magnitude of 8Ri 0 , such substitution implies that 8Rj 0 is always positive (Jekeli and 
Upadhyay, 1990, p. 10,984), while there is nothing dictating a fixed sign for 8Ri 0 , 
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despite the fact that Rj 0 itself is a non-negative quantity (see equation (2.143)). As it can 
be seen from equation (2.146), 8Rj 0 carries the same sign as cos|3, and the sign of cosP 
cannot be determined since it depends on the relative orientation of if 0 and the unknown 
vector 8i io. To the author's knowledge, this erroneous substitution appeared originally 
in Hajela’s derivations (1978, p. 5, equation 2.2), and has been adopted afterwards in a 
number of investigations including the recent study by Jekeli and Upadhyay (1990). The 
maximum magnitude of 8Rj 0 for a reference model complete to degree and order 8 is 
about 0.4 mgals for the STAGE mission considered by Jekeli and Upadhyay (ibid), 
hence 6R| 0 cannot be neglected in view of the 0.3 mgal total noise level that they have 
estimated for the (pseudo) observations Avi<> (ibid, p. 10,983). 

The existence of the centrifugal term 8Rj 0 in equation (2.142) makes the analysis 
of the SST observational system, in terms of the pseudo observations Av; 0 , practically 
untractable. In contrast, such a problem is not encountered if the primary GPS 
observables (carrier phase and pseudo-range) are used in a global dynamic solution for 
the simultaneous estimation of satellite orbital parameters and geopotential coefficients 
(Pavlis, E. et al., 1990). In the current application, which is of local character, to 
minimize the magnitude of 8Ri 0 , so that its omission can be justified, one may consider 
the following strategies: 

(a) Use of a higher degree reference geopotential model for the formation of Avj 0 . As 
it can be seen from Table 3, a model complete to Nmax = 20, implies a maximum 
expected value of 8Ri 0 of about 23 pgals, which may be considered negligible. The 
maximum degree of the reference model however, should be selected with due 
consideration to the ratio of the residual signal to the noise of the pseudo 
observations Avj 0 as was discussed before. 

(b) In a Multiple-High- Single-Low SST configuration, where as many as 8 GPS 
satellites are being tracked simultaneously by the low orbiter, one may select for the 
formation of Av; 0 pseudo observations, those GPS satellites for which the quantity 

fio = I £ fol COSY 

Pio 


(2.153) 
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is minimized (and/or does not exceed a pre-established threshold value). Such editing 
criteria should be considered in conjunction to other geometric requirements (such as 
minimum elevation angles of the GPS satellites). It is obvious that such an editing 
procedure introduces substantial complexity to the data processing algorithm. 

Finally, the noise Ej 0 of the pseudo observations Avi c is considered following the 
lines of Jekeli and Upadhyay (1990). As it can be seen from their analysis (ibid, p. 
10,976, Table 2), in the non-differential mode of observation, the major error source 
contributing to ejo arises from the frequency instability of the GPS satellite oscillator. 
This error source effectively cancels out if single differences are formed between the 
phase measurements to a GPS satellite as observed by the low orbiter and a ground 
receiver. In such case, the total error eu, of the residual acceleration along the line-of- 
sight, for the case of GP-B (10 cm orbit error, drag-free instrumentation), was estimated 
to be about 0.2 mgals. 

In the following error analysis however, it will be assumed that the SST 
configuration contributes information on the vertical component of the gravitational 
acceleration at altitude. The error of the vertical component of the acceleration may be 
approximated by multiplying Ej 0 by the vertical dilution of precision (VDOP) (Jorgensen 
1980) as described by Jekeli and Upadhyay (1990, p. 10,978). In their analysis it was 
found that for the Shuttle being the low orbiter, a representative value for VDOP was 
about 2. Adopting the same value for GP-B (which is a rather conservative assumption), 
one finally obtains 0.4 mgal error for the vertical component of the residual gravitational 
acceleration at GP-B altitude. 


CHAPTER III 


GLOBAL MEAN SQUARE ERROR ESTIMATION 

As it was discussed in section (1.3), the geopotential differences between points 
separated by ocean, whose geocentric positions are known with accuracy of a few 
centimeters, will be estimated by combining information from: 

(a) gravity disturbance measurements on the Earth's surface 

(b) gravitational acceleration "observations" at the altitude of a low Earth 
orbiter 

(c) a global geopotential model. 

The development of appropriate analytical tools that can be used to estimate the 
expected accuracy of the resulting geopotential differences, based on assumptions related 
to the spatial arrangement and the uncertainties of the input data, is the subject of this 
chapter. In practical applications, the geopotential differences will be estimated using 
least-squares collocation (Isc) (Moritz, 1980) since this technique can combine efficiently 
observations of different functionals of the gravitational potential, and assures that the 
resulting estimates are based on optimal use of the information contained in the 
observations. In addition, lsc can be used to assess the expected uncertainty of the 
estimated values, and thus offers one possibility for the means with which the current 
error analysis could by conducted. 

For given error properties of the input data, the (true) error of the estimated 
geopotential difference depends on the relative positions of the observation points with 
respect to the points to which the difference refers, as well as on the absolute positions of 
the latter. In that sense, the error of the estimated geopotential difference between points 
which lie on areas where the gravity field changes rapidly is likely to be larger than 
between points which lie on smoother areas of the field. The (squared) error estimate 
obtained from lsc on the other hand, represents a global average value, whose meaning is 
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as follows (Colombo, 1980). Keeping the relative positions of observation and 
estimation points fixed, subject the entire pattern to all possible rotations over the sphere. 
After each rotation measurements are taken, the lsc estimation is performed and the 
corresponding true estimation error is evaluated somehow. The global average of the 
squares of these errors, represents the squared error obtained from lsc. As a result of the 
use of homogeneous and isotropic covariances, if the relative geometry of observation 
and estimation points remains fixed, and so do the error properties of the data, the error 
estimate obtained through lsc will remain unchanged regardless of the absolute position 
and orientation of the whole observation/estimation pattern on the terrestrial sphere. 
Moreover, as it is well known, the lsc estimator, by its definition, ensures that this global 
mean square error is the minimum among the corresponding ones of any other linear 
estimator (Moritz, 1980, pp. 122-132). This property, as well as the ability of lsc to 
accommodate different data types and arbitrary spatial arrangements of 
measurement/estimation points, are responsible for the wide application of the lsc 
estimation technique. 

To benefit from the above properties however, one must take up the computational 
effort of forming and inverting covariance matrices whose dimensions equal the number 
of observation points. In actual implementation where the geometry of 
observation/estimation points is given and one is interested in the most rigorous 
evaluation of the estimates themselfs, as well as their expected errors, such computational 
effort is well justified. However, in error analysis studies the relative geometry of 
observation and estimation points is more or less a matter of assumption, and one is 
interested only in the expected errors of the estimates and not the estimated values 
themselves. In addition, if one strives for estimates which will not be significantly biased 
by uneven data distribution, it is necessary to impose some kind of requirements 
pertaining to the uniformity and spatial density of the measurements. Under such 
circumstances, if one is willing to accept additional assumptions that result in symmetric 
patterns for the geometry of the observation/estimation points, efficient techniques can be 
used to assess the expected errors of the estimates, resulting in large computational 
savings. Obviously, the reliability of the error estimates evaluated this way , depends on 
how well the assumed symmetric data arrangements compare with actual data 
configurations, at least in a global average sense. 
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The simplest way to assess the error of a geopotential difference, estimated on the 
basis of a global geopotential model and terrestrial- onlv measurements, is through error 
propagation based on truncation theory. As it will be seen however, this technique 
cannot be used if gravitational accelerations at altitude are to be included in the estimation. 
On the other hand, least-squares collocation using "ring averages" (Colombo, 1980, 
section 4.1) can accommodate both data types (either separately or in combination) and is 
only restricted by assumptions related to symmetries in the measurement pattern. In that 
sense, the use of ring averages provides a compromise which maintains the efficiency of 
integral formulas while incorporating the versatility of Isc. The analytical formulation of 
these techniques and their intercomparison are discussed next. In addition, the 
covariance models for the signals and the noise of the measurements which will be used 
in the numerical analysis are presented afterwards. 


3.1 Error Propagation Using Truncation Theory 

The use of truncation theory for the assessment of the global mean square error of 
geoidal undulations estimated from gravity disturbances measured over a cap centered at 
the computation point, and a global geopotential model, has been considered variously by 
Jekeli (1979) and Sjoberg and Fan (1986). In the following discussion the notation used 
by Despotakis (1987), for the corresponding case of cap integration of gravity anomalies 
(Stokes' kernel), will be adopted. 


Equation (2.102) of the previous chapter is written for convenience as follows: 


ifs, X ) -^|| H(\jr)D(S\ A/)dc 


(3.1) 


where: 

d(S', X') = [N b (8g + 8g A - e P + gi) - eJS', A/) . ( 3 . 2 ) 

The n* - degree and m* - order ellipsoidal harmonics of D and T are related by: 

D„(5', r)-(n+irrj5'.v) 


(3.3) 
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as it can be easily seen from equations (2.98) and (2.101). In addition, the solid angle 
corresponding to the cap centered at the computation point, is denoted by a c , so that: 


a ■ a c + (a - o c ) 


(3.4) 


and a kernel modification function Wj(\jr) is introduced, such that: 

H(\|f) = Hi(\|/> + Wi(\|f) ; 0<\|r^7t ( 35 ) 

and the modified kernel Hj(\|/) is defined through the last equation, given the definition of 
Wj(\|/). Due to (3.4) and (3.5), equation (3.1) becomes: 


a 

•iii « 


(y)D(5', X')do - 


4 n 


// - 

a-a c 

II w 


(V)D(5\ X')do 


(V)D(5 / , X')do+ JL 1 1 Wi(v)D(8', V)da 


(3.6) 


A function H; (y) is now introduced, defined by: 

Hi(V) - 

where y 0 is the semi-aperture of the spherical cap a c . Accordingly, equation (3.6) 
becomes: 


0 if 0 £ \\f <, \jT 0 

Hi(\|/) if % < v < 7t 



[H 1 (v)-H i (v)]D(5 , 1 X')da 
[h i (v) + W i (v)]D(5', X')dc . 


(3.8) 
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As long as the modification function Wj(\j/) is at least piecewise continuous in the interval 
0 £ v £ 7t (in which case Hi(y) and Hj(\|r) will also be piecewise continuous), it may be 

expanded in a series of Legendre polynomials, as: 


oo 

Wi(v) = X 2D T 1 W inPn(COS\K) 
n=0 Z 

(3.9) 

where: 


= I W i(\|/)P n (cosy) sinydy . 


K 

(3.10) 

Similarly, for Hi(y) and Hi(y) one has: 


oo 

H i(V)= X Qin(Vo)Pn(COSV) 

n=0 L 

(3.11) 

where: Qin(Vo)= 1 Hi(y)Pn(cos\|Osin\|rdY => 

Jo 


Qin(Vo) = I Hi(v)P n (cosv)sinvdy 



(3.12) 

and: 


oo 

Hi(y) =X X inPn(COSV) 

n=0 Z 

(3.13) 

X in = 1 Hi(v)P n (cosv)sinvd\ir 

Jo 

(3.14) 
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Making use of the decomposition formula of the Legendre polynomials (Heiskanen and 
Moritz, section 1-15), and the orthogonality of the surface harmonics, it can be shown 
easily that due to (3.9), (3.1 1) and (3.13), equation (3.8) becomes: 

©O oo 

t(S, X) = 4 £ [Xfa - Qin(Vo)] D„(8, X) + 1 X [Win + QM] D n (8, X) 

Z n=0 Z n=0 (3.15) 

where: 

D n (5, X) = X DJ5, X) . 

m=^n (3.16) 

Equation (3.15) is rigorously equivalent to (3.8) and provides the "frequency 
domain" counterpart of the latter. From (3.15) and (3.3) it can be seen that the quantities 
Xjn and should fulfill the relation: 

Xi n+ W m=^T ; n £0 (3.17) 

which provides the means of evaluating Xj n , without performing the integration (3.14), 
once Win has been evaluated. 

For the purpose of practical implementation, one has to consider that one part of 
the available information (the gravity disturbance measurements) represents "space 
domain" quantities, while the other part (the global geopotential model) is given in terms 
of spectral components. It is thus reasonable to seek a combination of equations (3.8) 
and (3.15) such that both kinds of information can be considered simultaneously in an 
efficient manner. Obviously, such a combination is meaningful only if the cap 
measurements are more detailed and/or accurate than what can be deduced for their values 
from the global model, and the cap does not extend to vjr = it. In that sense, one has from 
(3.8) and (3.15): 





H;(\|/)d( 5', X')da + 



[Win + QUVopnfS, X) 


(3.18) 
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whereby the integral term in (3.18) represents the cap contribution to T(8, X), implied by 
the modified kernel Hj(y), and the infinite sum represents the remote zone contribution to 
T(8, X) implied by the original kernel H(y) plus the cap contribution of the kernel 
modification Wj(y), i.e. 


^■StWin+cMvopis, x) = 

L n=0 


4)t 


Jj H(y)D(8',r)do+ 

a - ofc 


4lt 



(3.19) 


Equation (3.18) (which is rigorously equivalent to both equations (3.8) and 
(3.15)), forms the basis upon which a computational formula suitable for practical 
implementation can be developed. Obviously, different choices of the modification 
kernel Wi(y), yield different estimators of T(8, X), with varying properties. However, 
regardless of the choice of Wj(y), equation (3.18) states that to determine the mis value 
of T(8, X), requires continuous and errorless data D(8\ X') inside the cap of integration 
a c , and perfect knowledge of the spectrum D„(8, X) up to infinite degree. In practice, 
none of these requirements can ever be fulfilled; cap measurements can only be acquired 
at a finite number of discrete locations and are contaminated by observational noise, and 
in addition the knowledge of the spectrum extends to a finite degree and is imperfect. 
The errors that these imperfections induce to the estimated value T(S, X) of the disturbing 
potential are examined next, along similar lines to the derivations of Christodoulidis 
(1976). 

In practical implementation, the cap integration in (3.18), has to be replaced by a 
finite summation, since the function D(8', X') is not given in an analytic form but has 
only been sampled at discrete locations (Heiskanen and Moritz, 1967, section 2-24). To 
avoid biasing the result of this numerical integration, due to uneven distribution of the 
point measurements, one usually forms from the original point data, a set of area-mean 
values on a fine grid (e.g. 2’ x 2'), covering the entire cap ct c (Despotakis, 1987). It is 

assumed hereon that the cap data consist of area-mean values Dy, on an equiangular grid 
in terms of spherical distance y and azimuth a (Ay = Aa), centered at the computation 
point P(8, X). Note that y and a are evaluated from ellipsoidal coordinates (8, X) and 
(S', X'). The indices (i, j) identify the location of the compartment to which the area- 
mean value refers, in a two-dimensional array where: 


i = 0, 1,2, Nr - 1 


; j = 0, 1, 2, .... 2N - 1 


71 

(3.20) 


As illustrated in Figure 9, the number of "rows" (or "rings") of area-mean values around 
the computation point is Nr, while each ring contains 2N compartments, where: 

N = — — . (3.21) 

Aot 


j = o 



Figure 9. Arrangment of Cap Data Around the Computation Point P. 

In addition, the estimated area-mean value D-j, differs from the corresponding true value 
Djj, by the true error eD^ , so that: 

+ • (3.22) 

As far as the available spectral information is concerned, a finite number of harmonics, 
D^(8, X), up to maximum degree M, is assumed to be given. These are contaminated by 
commission error eDS(8, X), so that: 
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(3.23) 


DnfS.X^DfexJ + eD^S.X) ; n = 0, 1, .... M 

where Dj(8, X) represents the true value of these harmonics. Based on and D^(8, X), 
the following estimate of the disturbing potential is evaluated: 

Nr-l 2N-1 [ f M 

f( 8 , X)= 2 L X X DC ij| I Hk(v)da + [Wk„ + QUVo)]Dn( 5 , X) (3.24) 
471 i=0 j=0 J) z n=0 

where ay is the solid angle (area on the unit sphere) corresponding to the (i, j) th 
compartment, and the subscript k was used above to discriminate different kernel 
selections. The true error of the estimate t( 8, X) is given by: 

elfs, X)= t(8, X) - t(8, X) , (3.25) 


so that, taking into account (3.18), and (3.22) through (3.25), one has: 

et(8, X) = ei(8, X) + e 2 (s, X) + e 3 (s, X) + e 4 (s, x) ( 3 . 26 ) 

where: 


ft Nr-l 2N-1_ t [ 

ei(8.X) 1 1 H k ( V p(5', Xjdo - ^ X i D ijJ I Hjy) 

o c °y 

, . Nr-l 2N-1 

£ EDiJJ Hk(¥) 

Oij 


da 


da 


M . 

£3(5, x) = 4 X [Wkn + QidM eE^(5, X) 


n=0 


e 4 (8, X) = ^ £ [Wkn + Qkn(Vo)]D n (8, X) 

2 n=M+l 



(b) 




(3.27) 


(c) 


(d) 


I 
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The following nomenclature will be used for the above error components 
(Christodoulidis, 1976): 

El : sampling (or discretion) error of the cap integration 

£2 : propagated error of the cap data 

£3 : commission error of the geopotential model 

£ 4 : omission error of the geopotential model. 

It should be emphasized that eDij refers to the area-mean value of the (i, j) th 
compartment, and not to the original point measurement, although in practice, for 
compartments as small as 2 ' x 2 ', the estimated error of such area-mean value differs 
very little from the corresponding error of a point measurement. 

The error component £i( 5, X) is considered next in detail. The following notation 
is introduced: 

£ 1 ( 8 , A.) = a( 8 , X) - b( 8 , A.) (3.28) 

where: 


»(S,X)=^j| H i t ( vp(S',V)do 


Oc 


Nr-1 2N-1 [ f 

b ( 8, X) =£ 2 2 D ijJJ H k(v) 

J " °ij 


(a) 


da (b) 


From equation (3.15) one has immediately: 


a(5, X) —4 ^ [Xfcn - Qkn(Vo)J^n(8, X) 

L n=0 


(3.29) 


(3.30) 


while for b( 8 , X), employing in a discretized manner the same technique used before to 
derive (3.8) from (3.6), one can write: 


b(8, X) = bi(5, X) - b 2 (8, X) 


74 

(3.31) 


where: 


N-l 2N-1 ( [ 

bi(S.X)=£I iDuljHiM 

J " Oij 
N-l 2N-1 ( (_ 

b 2 (8A)= s L£ SDijljHkM 

J " Oij 


da (a) 


(3.32) 


da (b) 


Consider the term bi(8, X) first; due to equation (3.13) and the decomposition formula: 
P„(cos V )= Yj8,x)Yj8',r) . (3 33) 

equation (3.32a) becomes: 


r r 


N-l 2N-1_ 

b,(8, X) 4*? ^° ij 


i=0 j=0 


£2nil yJs, A.) yJs', A.') do 

L^O 2 2n+ 1 m =-n J 


Oij 


N-l 2N-1 o. n . [ [ _ . . 

bi(8, X) = ^>j S ^ kn S Y nm (8, X) 1 1 Ynm(s , X ) 

8JC i=0 j=0 n=0 m= n J J 

Oij 


da 


(3.34) 


The coordinates of the variable point of integration (S’, X’), are related to the integration 
variables spherical distance y and azimuth a with respect to the computation point 
P(8, X), by the well known formulas (Heiskanen and Moritz, 1967, p. 113): 


cosy = cosScosS' + sin8sin8'cos(X' - X) 


sin8'sin(X - X) 

tana= ; 

sinScosS' - cos8sin8'cos(X' - X) 



(3.35) 


and the integral over Ojj will be denoted by IY nm (i j) so that: 
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t r /Vi+Ay i-aj+Aa 

1 1 Ynm (S', X')da = I I Y nm (8', X') sinydyda 
a ij ^ 


- IYnm(l,j) • (3 3 6 ) 


Re-arranging the order of the summations in (3.34) one has, due to (3.36): 


bi(5,X) = l f;x kn £ Y nm (8,X) 


n=0 m =-n 


N-l 2N-1 _ 

4. Z Z DijIYnm(i» j) 
. i=0 j=0 


(3.37) 


The expression inside the brackets in (3.37) is readily recognized to represent a 
quadrature formula (Colombo, 1981a), according to which an estimate E^ m of the true 
spectrum D nm of D(8, X) may be obtained, given a global set (N x 2N) of errorless area 
mean values Djj, on the equiangular grid A\|/ = Aa with pole the computation point 
P(5, X). Hence: 

N-l 2N-1 _ __ 

6 nm= Z DijJYnmaj) . 

i=0 i=° (3.38) 

It should be mentioned here that the coefficients E|, m obtained from (3.38) refer to 
the coordinate system associated with the (8, X) grid, although the area-mean values Dy 
are given on the (\j/, a) grid with pole the computation point P(8, X). Although in theory 
the two grids do not have to coincide, this incompatibility makes the practical evaluation 
of the quadrature formula (3.38) extremely inefficient, basically because the computation 
of integrated Legendre functions (upon which the evaluation of IY^fi, j) depends), 
cannot be accomplished using efficient recursions as those derived by Paul (1978). 
However, equation (3.38) is used here as an intermediate step to provide insight and aid 
in the derivation of the final expression for the sampling error. As it will be seen, the 
final formula for the sampling error does not require for its evaluation the actual 
implementation of (3.38). Due to (3.38) equation (3.37) becomes: 


bi(5, X)= X kn D n (8, X) 

L n=0 


(3.39) 


where: 
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D„( 8 ,A) = Z D nm Y nm ( 8 , A) . 

m=-n (3.40) 

Following the same procedure for b 2 ( 8 , A.), and taking into account equation (3.31), one 
obtains for b(5, A): 

oo 

b(6, X) = 1 X [X ta - Qkn (Vo)] D n (8, X) . 

z n=0 (3.41) 

Hence, due to (3.28), (3.30) and (3.41) the sampling error ei(8, A.) finally becomes: 

oo 

£l(s, A) = \ X [Xkn - Qkn(Vo)]Sn(5, X) 

1 n=0 (3.42) 

where: 

S n (5, A) = D n (8, X) - D n (5, X) = I [Dm, - D nm ] Y^S, X) 

m=^n (3.43) 

and S n (8, X) represents the n* - degree surface harmonic of the sampling error associated 
with the quadrature formula (3.38). 

Proceeding along similar lines, one has for the propagated error £ 2 ( 8 , X): 

~ n _ N-l 2N-1 _ _ 

e2(S,X) = lX[ X kn-Qkn(Vo)]Z Y nm(8,X)^-X Z eD ij ^(i, j) . (3.44) 

n=0 m=-n |_ 47C i=o j=0 

One may view the true errors eDjj of the area-mean values as discrete samples of an 
error function defined over the full solid angle (unit sphere). To be more precise, the 
totality of the eDjj values represents one realization of a stochastic process on the unit 

sphere (Moritz, 1980, p. 279). This error function can be expanded in surface 
(ellipsoidal) harmonics, and the coefficients E nm of such expansion may be approximated 

by the following quadrature expression, given the discrete samples eDjj : 

N-l 2N-1 _ __ 

£nm ~ ~ Z. Z lYmnO, j) 

4 ^Pn i=0 j=0 


(3.45) 
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where p n is the Pellinen smoothing factor of degree n, which corresponds to a spherical 
cap having area equal to the area of an equiangular compartment in the ring with i = N/2 
(the "equator" of the (\\f, a) grid). It is acknowledged here that the use of smoothing 
factors p n which are independent of \|/ (or equivalently of the ring index i) is only an 
approximation, since the area of the equiangular compartments on the (\y, a) grid varies 
with \jr as it can be seen from Figure 9 (detailed discussion of this aspect is given in 
Katsambalos (1979)). Using a complete (N x 2N) set of area-mean gravity anomalies on 
an equiangular grid, it has been verified numerically that the spectrum implied by a 
quadrature formula of the type (3.45) (ring-independent p n ) differs from the spectrum 
implied by a quadrature formula using ring-dependent p n factors, by about 15% near the 
degree corresponding to the Nyquist frequency (n = N). In view of the fact that the error 
properties assigned to the data are to a significant extent a matter of assumption, the 
approximation introduced in (3.45) by using ring-independent p n factors appears to be 
acceptable. Accordingly, equation (3.44) becomes: 

oo 

£2(8, X) - 1 X Pn[ X kn * Qkn(Vo)] £n(5, X) 

2 n=0 (3.46) 


where: 


£n(6,X)=X enmYnmM) 

m=-n (3.47) 

The formulation suggested here to model the sampling (or discretion) and the 
propagated errors is slightly different than the one originally proposed by Christodoulidis 
(1976) and adopted by Despotakis (1987). Christodoulidis (1976) examined alternative 
techniques to model the sampling error, and concluded that a computationally manageable 
and accurate enough model would be (in the current notation): 

oo 

£i(S, X) = ± X [X kn - Q kn (Vo)]D n (5, X) . 

2 n=N+l (3.48) 

(The fact that the analyses cited above were made for geoidal undulations obtained from 
integration of gravity anomalies is immaterial, since the four error sources identified in 
equation (3.27) are present in both integral formula applications.) The model (3.48) is 
based on the assumption that the data inside the cap of integration contain spectral 
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information only up to degree N, which is related to the size of the data compartment by 
the well known rule of thumb given in (3.21). This is one possible approximation; 
however, equation (3.48), rigorously interpreted, represents the sampling error for the 
case of continuous coverage inside the cap with band-limited data (containing spectral 
information only up to degree N). Such an approximation of the real-world situation 
(where continuous coverage can never be achieved, and the spectral content of the area- 
mean values is rather difficult to assess), is less realistic, and does not yield a 
significantly simpler model, than the model (3.42) proposed here. 

As far as the propagated error is concerned, the Pellinen smoothing factors have 
been introduced here to account for the fact that area-mean values are used for the cap 
integration, while error covariance models in practice usually describe error properties of 
point data. Despotakis (1987, equation 2.26) altered the original formula for the 
propagated error as given by Christodoulidis (1976, equation 150), by truncating the 
error spectrum at degree N, effectively assuming that all error contribution above this 
degree is smoothed out from the area-mean values. In the model (3.46) proposed here, 
the error tapers off gradually as the degree n increases, since P n * 0 as n — > 

According to the above, the expected global rms value of each error component 
can now be derived. The following definitions are introduced first for notational brevity: 

XQkn = Xfcn “ Qkn(Yo) 00 | 

(3.49) 

WQkn = Wkn + CWVo) (b)J 

where it is understood that both XQkn and WQkn depend on the semi-aperture \|/o of the 
cap. Taking the total average I (Moritz, 1980, p. 100) of the square of each error 
component et(i = 1, 2, 3, 4), one has due to the orthogonality of the surface harmonics: 

oo 

rms^e,) = E[e^8, X)Ml XQU 

4 n =o (3.50a) 


and similarly: 
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rms 2 ^) = \ X (pnXQ k „) 2 a n (b) 


n=0 


M 


rms2(e3) = jX WQknEdn 


(c) } 


n=0 


(3.50) 


rms 2 {e4) = 2 X w QL d n (d) 

4 n=M+l 

where the linearity of the operator E was used, and the following notation was 
established: 


s n = e[s 2(5, X)]=M[sfe X)] (a) 


a n = E[46, x)] 

ed n = E[eD^(8, x)] 


(b) 

(c) 


(3.51) 


d n = H [d 2(5, X)] = M[D2{5, X)] (d )/ 


The homogeneous and isotropic space averaging operator M is defined by (Moritz, 1980, 

p. 82): 


m in 

t=0 


(°)sin8d5dXda 


(3.52) 


In addition, the total average of the products e*(A).e*(B) (£ = 1, 2, 3, 4) of each error 
component for two points A and B separated by spherical distance \jr<i is given by 
(Christodoulidis, 1976): 


rms 2 [ei(A).ei(B)] = E[ei(8 A , X A ).ei(5 B , X B )] 


= \ X x QL s n p n(cos\ir d ) 

4 n=0 


(3.53a) 


and similarly: 
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rms 2 [e 2 (A).e 2 (B)] = ^ X (PnXQtnf o n P n (cosy d ) (b) 

n=0 

M 

nn S 2[e3(A).e3(B)] = j X W^ed^cosyi) (c) 

4 n=0 

t 

oo 

rms 2 [e4(A).e4(B)] = 7 X w QLidnPn(cosy d ) (d) 1 

4 n=M+l 


0 . 53 ) 


Since the error components identified in (3.27) originate from independent causes, 
it is reasonable to assume that these errors are uncorrelated (it needs to be assumed here 
that the global geopotential model has been derived independently of the cap data). In 
such case, the expected global rms error of the geopotential estimated by (3.24), is given 
by: 


rms(ef) = ^ 


rms 2 ^) 


I1/2 


(3.54) 


while the expected global rms error of the geopotential difference cTab, between the 
points A, B separated by spherical distance yd is given by (Christodoulidis, 1976, p. 43): 


rmsfeTAB) ■ V7 ( X rms^e/fA)] - X ms2 fo( A )- e /( B )]| 

u=i / =1 I 


(3.55) 


The above formulation enables one to estimate the expected global rms errors in 
geopotential and geopotential differences, obtained from the estimator (3.24), once 
appropriate models have been established for the degree variances s n , On. £dn and d n . 

As far as the selection of the kernel modification function is concerned, two 
alternative principles may be followed for the definition of Wjt(y): 


(a) Deterministic approach: The selection of Wk(y) here is made in such a way, that the 
resulting kernel function that corresponds to the remote zone contribution ("truncation 
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kernel" - see equation (3.19)), possesses eigenvalues which converge to zero more 
rapidly than those corresponding to the unmodified kernel H(\|/) (Jekeli, 1980). In this 
manner, one attempts to reduce the error arising from the lack of detailed data outside the 
cap of integration, taking also advantage of the information provided by a global 
geopotential model. The increased convergence rate of the eigenvalues of the truncation 
kernel is accomplished analytically, by removing die discontinuity of this kernel (Meissl's 
method) or the discontinuities of the kernel and its derivatives (Molodenskii's method) at 
\|r = Y 0 , as discussed in detail by Jekeli (ibid). Such modifications are made without 
considering the error properties of the cap data or of the available geopotential model. 

(b) Stochastic approach: The determination of Wk(\|/) here is accomplished numerically, 

/s. 

by imposing the condition that the resulting global rms error of T is minimum. Such 
condition yields a linear system of equations for the eigenvalues Wkn of the modification 
kernel. This approach has been put forward by Colombo (1977), who considered the 
minimization of the truncation error only, and was developed further by Sjoberg (1986) 
to account for all error sources identified in equation (3.27). The eigenvalues Wkn 
determined in this manner depend on the assumed error properties of both the cap data 
and the geopotential model. 

Despotakis (1987) intercompared Stokes' kernel and its modifications according to 
Meissl's, Molodenskii's and Sjoberg's techniques, in a global rms error analysis fashion, 
as well as in actual geoidal undulation computations. His analysis indicated that in a 
global average sense and for cap sizes smaller than about 5°, the computationally simpler 
technique of Meissl is expected to be almost as accurate as the more demanding 
techniques of Molodensky and Sjoberg. In Despotakis' actual computations, Meissl's 
and Sjoberg's techniques yielded results of practically the same quality (ibid, p. 97, Table 
28). 

For the following analysis the original kernel of Hotine (k = 1) and its 
modification according to Meissl's suggestion (k = 2) will only be considered. For these 
cases one has: 


Unmodified Hotine's kernel (k = 1) 

Wi(y) =0 ; Hi(\j/) = H(\jr) 0<y^Jt 


(3.56) 
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and the eigenvalues of the kernels Wi(y), Hi(y), and Hi(y) are given by: 


Wj n = 0 


Xln = n^T 


(a) 


(b) 


Qin(Vo) = Qn(Vo) = I H(y)P n (cosy)sinydy (c) 

JVo I 


(3.57) 


for n £ 0. The coefficients Q n (y 0 ) should not be confused with the corresponding 
Molodensky truncation coefficients referring to Stokes' kernel. 

Meissl's modification (k = 2) 


W2(y) = H(y 0 ) - Hq ; H2(y) = H(y) - Hq 0<y <7t 
and the corresponding eigenvalues are: 

2Ho if n = 0 


(3.58) 


!n= i 


W2n= I H 0 P n (cosy)sinydy => W2n = 


0 if n>0 


(3.59a) 


X 2n = 


n + 1 

_ 2 _ 

n + 1 


-2Ho if n = 0 


if n>0 


1 


Q 2 n(v|/ 0 ) = I [H(y) - Ho] P n (cosy)sinydy = QAy 0 ) 


P ” 


(3.59b) 
(cos\|/)sinydy . 


The last equation, due to the recurrence relations of the Legendre polynomials (see 
Appendix A, equation A. 15), can easily be reduced to: 


{ Ho(l + cosy 0 ) if n=0 

- 

[Pn+i(cos\|/o) - Pn-i(cosyo)] if n >0 


(3.59c) 
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The above expressions for W2 n , ^2n> and Q2 n (Vo) are strictly valid for 0 < \|f 0 ^ tc, due 
to the singularity of Hotine's kernel at \jr = 0. The case \|r 0 = 0 implies that no cap data 
are used in the determination of T, hence only the commission and omission error of the 
geopotential model contribute to the error of T. In such case the sum of the eigenvalues 
Wkn plus Qk n (\|/ 0 ) (which appear in equations 3.27c and 3.27d), always equals 2/(n + 1), 
regardless of the selection of the kernel modification W k (\|/). 

From the numerical point of view, it is obvious from the previous formulation, 
that only the evaluation of Q n (Vo) is a computationally complicated process. Jekeli 
(1979) developed an efficient recurrence relation for the evaluation of Qn(Vo) (Vo * 0). 
His technique was expanded in this study and a similar recurrence relation was developed 
for the more general case of Pizzetti's extension of Hotine's kernel, i.e. for the kernel 
H(R/ r , \|0 where r > R. The detailed derivation of the recurrence relation for the 
truncation coefficients Q n (R/ r , \|/ 0 ), corresponding to this kernel, is given in Appendix A. 
In case one postulates the absence of the zeroeth- and first-degree harmonics from the 
gravity disturbance, the kernel H(R/ r , \|/) requires corresponding modification (Jekeli, 
ibid), and the resulting kernel H*(R/ r , y) (i.e., H(R/ r , \|/) less its zeroeth- and first- 
degree harmonics), implies a different set of truncation coefficients Qn(R/r> Vo)- The 
derivation of a recurrence relation for Qn(fyr> Vo) is given in Appendix B. 

The zeroeth-degree coefficients Qko(Vi)» i = 0, 1, 2, ..., N r provide an efficient 
way for the evaluation of the discretized integral over the cap, as long as the cap 
measurements are defined on the (\|r, a) grid. Denoting the cap contribution to the 
disturbing potential by Cap(\j/ 0 )» so that (see equation 3.24): 

Nr-l 2N-1 _ 

Cap(Vo)=^rX X D u 

1 i=o j=o 

(3.60) 

one can easily see that due to (3.12): 

H k (\|/)da = Aa[Q ko (Vi) - Qko(Vi-n)] 

Oij 

so that (with obvious notation for Qko)- 




(3.61) 
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and this equation demonstrates the fact that the cap contribution to the disturbing potential 
is the weighted sum of the ring-average "observations" Df, with weights the quantities 
AQko(i)- Equation (3.64) represents a rigorous evaluation of (3.60), and is a 
consequence of the isotropy of the kernel Hk(y) (no dependence on azimuth) and of the 
particular selection of the data grid (y, a). The form of equation (3.64) is valid even if 
the "width" of the rings (i.e. yi+i - yO and/or the number of compartments per ring (i.e. 
Aa) vaiy with the distance y from the computation point, as described in Heiskanen and 
Moritz (1967, p. 120), provided of course that Di and AQko(i) are evaluated 
accordingly. As it is well known (ibid, p. 120) the disadvantage of using the (y, a) grid 
to register the data, is the need to re-evaluate Djj (or Di ) as one moves from one 
computation point to another. The alternative of course, is to register the data with 
respect to the (8, X) grid (ibid, p. 117), but then the integration of the kernel function 
over the data compartment has to be performed numerically. Provided the original point 
measurements are available (so that one can decide upon which grid to use to register the 
mean values) and if the computation points are few in number and randomly distributed, 
the (y, a) grid is inherently better suited for cap integration. 
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3.2 Error Assessment for the Over-Determined Boundary Value Problem 

The formulation given in the previous section for the error estimation of the 
disturbing potential (or the disturbing potential difference), obtained from cap integration 
of gravity disturbances and a global geopotential model, is applicable only in case the 
computation points lie on or outside the sphere to which the gravity disturbance data 
refer. This is due to the fact that Hotine's integral formula (upon which the entire 
derivation was based) is the solution of Neumann’s boundary value problem for the 
exterior space of a sphere. The estimation of the error of disturbing potential on the 
Earth's surface, obtained from gravity disturbance measurements at altitude requires 
alternative treatment. In addition, if observations at altitude are to be used in combination 
with terrestrial measurements, then the boundary value problem in question becomes 
over-determined, as opposed to the uniquely-determined Neumann's problem solved by 
Hotine's integral. The statement of such a problem, for the case where the known 
boundary surfaces are concentric spheres may be given as follows (see also Figure 10 for 
notation definitions): 

"Determine a function T, harmonic in the infinite region outside the sphere (O, Rt) and 
regular at infinity, if its normal (radial) derivatives (-3T/9r) obtain prescribed values on 
the two concentric spheres (O, Rt) and (O, Rs) where Rs = Rj + h and h > 0." 

Since the boundary values on (O, Rt) alone, are sufficient to determine T outside 
(O, Rt) uniquely (through the solution of Neumann's problem), and since the same 
holds true for the boundary values on (O, Rs) and the space outside (O, Rs), existence of 
a unique solution to the over-determined problem requires a compatibility condition to be 
fulfilled between the two sets of boundary data. This condition can easily be obtained 
from the following considerations. Let Dj and D§ denote the boundary data on the 
spheres (0,Rt) and (0,Rs) respectively. Let also P(Rs, 0, X) be an arbitrary point on 
the surface of the sphere (O.Rs). Then, due to equation (2.67) one has: 


Tp = ff JJ H(y) Ds(0', V)d<J 


(3.65) 
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while due to equation (2.64): 


Tp = 


Rj 

4tc 



D^0', X')da 


(3.66) 


Hence, for a unique value of Tp to exist, Dj and D§ should fulfil the compatibility 
condition: 


II [h(|i , v) r t d^(8', V) - HM RsD|(e'. V) 


da = 0 


(3.67) 


for every point P(Rs. 0. X) on the sphere (O, Rs). In that case, due to Stokes' theorem 
(Heiskanen and Moritz, 1967, pp. 17-18), T is also uniquely determined outside the 
surface of the sphere (O, Rs). Due to the relation (2.73) between the radial derivative of 
Hotine's kernel and the kernel of Poisson, the condition (3.67) may also be replaced by 
the equivalent set of conditions: 


R S D|(R S , 0, X) = RT t R *> ™ - - - da 


r st L 

s 4ji 


ll D|(6', X' )da = r t^ j"j" D ^(0 / > X')da 


(a) 


(3.68) 


(b) 


Note that (3.68a) by itself is equivalent to a condition similar to (3.67), but with its right- 
hand side equal to a constant. Condition (3.68b) enforces such constant to be zero. The 
condition (3.67) (or equivalently the set of conditions (3.68a, b)) would have been 
fulfilled as an identity, if the measurements Dj and D$ were free of errors. In the 
presence of observational noise however, to enforce a unique solution to the over- 
determined bvp, one could use either (3.67) or (3.68a, b) and set up a least-squares 
adjustment of the form F(La) = 0 (condition equations), for Dj and D§. The result of this 
adjustment would be a set of Dr and Ds which would fulfil the compatibility condition 
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and would be accompanied by a-posteriori estimates of their error properties. The 
a-posteriori errors of E>r could then be used in the truncation theory formulation of the 
previous section, to estimate the error of geopotential (T) or geopotential difference (AT). 
In this manner the additional information provided by the measurements at altitude, for 
the estimation of T or AT on the Earth's surface, can be taken into account. However, 
the above procedure requires the measurements Dr and Ds to provide global coverage, so 
that (a discretized form of) (3.67) or (3.68) can be used to set up the condition equations. 
The problem at hand however, is of local nature. Terrestrial measurements are given in 
the caps centered at the computation points A and B, while the measurements at altitude 
are given in the caps centered at Sa and Sb as illustrated in Figure 10. 


Mean sphere of satellite 
observations (r = R$) 



Figure 10. Geometry of the Measurement and Estimation Points. 
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One way to circumvent the inability to use integral formulas for downward 
continuation, is to perform the error analysis in two steps: 

(a) Using the altitude generalized truncation coefficients derived in Appendix A, one 
can estimate the error eAT s J Sb of ATs a s b when the latter is estimated from 
terrestrial measurements via upward continuation (Pizzetti's extension of Hotine's 

c 

integral). This error may then be compared to the corresponding error eAT SaSb 
obtained when ATs a s b is estimated from the measurements at altitude. 

(b) Significant contribution of the observations at altitude to the accuracy of ATab is 
provided when eAT s J Sb > eAT| A s B . To quantify this contribution one may use 
least-squares collocation where the signal to be predicted is ATab and the two 
independent input signals are: ATs a s b obtained from data at altitude and ATab 
obtained from terrestrial measurements, each input signal accompanied by the 
corresponding error as obtained from truncation theory. 

Obviously the motivation for using the above procedure is to take advantage of the 
simplicity of error propagation through truncation theory on one hand, and of the "built- 
in" ability of least-squares collocation to perform the downward continuation (Moritz, 
1980, p. 97) on the other. However, as it will be shown next, these properties may be 
exploited in a more efficient way, if least-squares collocation with "ring-averages" is used 
for the error analysis. In such case all measurements are used simultaneously in the error 
analysis and the error estimate eATab is obtained in one step. 

3.3 Least-Squares Collocation Using Ring Averages 

Consider the observations in the caps centered at A, Sa (see Figure 10) arranged 
in the vector 4 a and similarly those in the caps centered at B, Sb in the vector dB- These 
observations consist of a signal part denoted z and a noise part denoted n, so that 
(Colombo, 1980): 
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The signal to be estimated from the data Ua and dB is the geopotential difference ATab. A 
linear estimate of this signal has the general form: 

ATab = fd (3.70) 


where: 



and the estimator f can always be considered as composed of two parts: 



(3.71) 


(3.72) 


Hence, the estimate ATab becomes: 

ATab= fBdB-fldA (3.73) 

and its error is given by: 

eAT ab = (T b - T a ) - (f5 dB - fllU) • (3.74) 

Applying the total averaging operator E (Moritz, 1980, p. 100) to the square of the 
estimation error one has: 

^(eAT ab) 2 ] = S(T B - T a J 2 + fj dsdjfe + - 2(T B - T A )dJfB 

+ 2(T B -T A )dIfA-2fidAdIfe] • (3-75) 


Assuming no correlation between signals and noise, i.e. 
HTn T ) = B(zn T ) = 0 


(3.76) 


and denoting: 
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M(Ti') = Ctt 

(a) 

MCIz 1 ) = C Tz 

(b) 

M(zz T ) = Czz 

(c) 

Etan 7 ) = D 

(d) 


(3.77) 


one has from equation (3.75): 

^(eATabF ] = 2Ctt (A, A) - 2Ctt (A, B) + fjjtCz^B, B) + D(B, B)]£b 
+ fftCWA, A) + D(A, A)]£a - 2[Ctz(B, B) - C Tz (A, B)]fe 
+ 2 [C Tz {B, A) - C Tz (A, A)]£a - 2fJ[C zz (A, B) + D(A, B)]fe (3 78) 

From all possible choices of linear estimators f one seeks now such an f that minimizes 
^(eATabF ]• Imposing the conditions: 

(3?9) 

one arrives, due to (3.78), to the following linear system of equations in and £b: 
$C zz (A, A) + D(A, A)] - fSCzzfA, B) + D(A, B)] + C T2 (B, A) - C T2 (A, A) = 0 
- fitCd(A, B) + D(A, B)] + fSC^B, B) + D(B, B)] - C Tz (B, B) + C T2 (A, B) = 0 

(3.80) 

The following assumptions are now made: the observations sIa and £ta have the 
same configuration and are characterized by the same error properties. In addition, no 
correlation exists between the noise ha and hb- If one seeks an estimate of ATab which 
should not be significantly biased towards either endpoint of the baseline AB, it is 
reasonable to require that at both endpoints A and B, the available measurements have 
similar configuration and error properties. Hence, the above assumptions are not 
unrealistic from the practical point of view. Under these assumptions one has: 


C Z2 (A,A) = C ZZ (B, B) 

(a)] 

Ctz(A, A) = Ctz(B, b) 

(b) 

C T z(A, B) = Ctz(B,A) 

(c) 

D(A, A) = D(B, B) 

(d) 

D(A, B) = 0 

(e)/ 


and the solution of the system (3.80) can easily be found to be: 

£a = £b (3.82) 

£a = [C ZZ (A, A) - Czz(A, B) +D(A, A)]*1[C£ (A, A) - C? x (A, B)] . (3.83) 

Substituting (3.81) - (3.83) in (3.78) one finds that under the previous assumptions the 
global mean square error of ATab is given by: 

^(eAT abF] = iCrr(A, A) - Ctt(A, B) - U T C 4 u] (3.84) 

where the vector U and the matrix C are given by: 

U = CxzfA, A) - Cf^A, B) 

C = CzzfA, A) - CzzfA, B) + D(A, A) 

Although the symmetry assumptions (3.81) significantly simplify the computation 
of the global mean square error of ATab. the formation of the covariance matrix C in 
(3.85b) requires computational effort proportional to the square of its dimension which is 
equal to the number of observations. The formation of C, as well as its inversion, can be 
a very demanding computational process. On the other hand, as it was shown in section 
3.1, the isotropy of Hotine's kernel in conjunction with the suitable selection of the (y, 
a) grid to register the cap data, resulted in very efficient formulas for the computation of 
ATab as well as its error, in the case of the uniquely determined bvp and the integral 
formula approach. In that case, as it was shown in equation (3.64), the cap contribution 
to ATab was given as a weighted sum of ring-average "observations". It is thus 
appropriate to consider the simplifications of equations (3.85a, b) if instead of the 


(a) | 

(b) | 


(3.85) 


92 

original measurements, ring-averages are to be used in the least-squares collocation 
estimator. Obviously, in such case the number of "measurements" reduces drastically 
and this brings significant savings in the formation of the covariance matrices (Colombo, 
1980). To follow this approach the analytical expressions for the covariances between 
ring-average "measurements" have to be developed first. This is done next 



Figure 1 1. Geometry Associated with Ring- Averages. 

Let t and s be two gravimetric quantities. The homogeneous and isotropic 
covariance between t and s (both considered as point values) may be expressed as: 


COV(t P , Sr) = M(tp Sr) = X °ts( n ) P n(cOSypR) 

n=o 


(3.86) 


where the averaging operator M was defined in (3.52) and <Jts( n ) (n = 0,1,2...) is the 
power cross-spectrum of t and s (i.e., Ots( n ) is the global average under M of the 
products of the surface harmonics t n (5, X) times s n (5, X)). The ring-average of the 
quantity s is defined as (Jekeli, 1989): 


s(Vqr) = 


.-Lf 

271 

^0Cqr=O 


S(VQR, CtQR)d(XQR 


( 3 . 87 ) 
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where the notation can be read from Figure 11. The covariance between the point value 
tp and the ring-average value s(\j/qr) is given by: 


covftp, s(\|/q R )] =M 




O 


2n 

s(Vqr, Oqr) dotQR 


2n 


I 


In 

M[tpS(YQR, (XQr)] dotQR 


(3.88) 


Thus the required covariance is the average over the azimuth <Xqr of the covariance 
between tp and Sr, where the point R now moves on a circle of semi-aperature \|Tqr, 
centered at Q. Substituting (3.86) in (3.88) one has: 


covft P> s(\|/Q R )] 



P„ (cos ypr) dotQR 


(3.89) 


where \|/pr depends on \j/pq, \|/qr and ocqp (which are constant), as well as on the variable 
azimuth ccqr. From the decomposition formula (Heiskanen and Moritz, 1967, equation 
1-82) one has: 

P n (cOS\|Tp R ) = P n (cOS\J/iQp)P n (cOS\J/qr) + 

+ 2 X [Rnm(VQP. OtqpjR^VQR. Oqr) + 

+ Snm(VQP» ttQPjSnmfVQR. «qr)] (3.90) 


where R n m and S n m are defined in (ibid, equation 1-67). The integral in (3.89) can be 
written as: 



P n (cosy PR ) cos(0.aQR) d(XQ R 


(3.91) 


so that upon replacing in (3.91) the expression (3.90) for P n (cos\yp R ), due to the 
orthogonality of sines and cosines in the interval [0, 2 jc], one has: 
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O 


P n (cosVPR) dctQR = 2n P^cosVpq) Pit(cos\j/qr) 


and (3.89) finally becomes: 


(3.92) 


cov[t P , s(\|/qr)] = X <y ts (n)Pn(cosVQR)P n (cosVpQ) . 

n=o (3.93) 

Following similar lines, the covariance between two ring-average values t(\j/py) (where 
\|/pv is the semi-aperture of a ring centered at P) and s(Vqr) can easily be found to be: 


cov[t(v PV ), s(v QR )] = X a ts(n)P„(cosvpv)Pn(cos\|fQR)P n (cosvpQ) . 

n=o (3.94) 

The general expressions (3.93) and (3.94) can now be specialized for the current 
application. If d n denotes the degree variance of the point value of gravity disturbance 
(more precisely of -dfT/dr), referring to the surface of a sphere of radius R = 6371 km and 
8g denotes ring-average value of gravity disturbance, then (see also section 3.4): 


oo 

Ctyi{P, Q) = X (t^yF (a) 

oo 

Cr,5 t (P, Q R ) = X ( n ^ ^ j" Pn(cosv QR )P ii(cosypq) (b) 

t ,5,(P V - QR) = X dn|^p 2 Pn(cOSVpv)Pn(cOS\|rQR)P n (cOSVpQ) (c) 


(3.95) 


Note that the summations in (3.95) start from n = 2, i.e., it is assumed that the disturbing 
potential does not include zeroeth- and first-degree harmonics. 

It should be emphasized here that equations (3.93) through (3.95) refer to rings 
formed by point measurements and not to zones of area-mean values in equiangular 
compartments as it was the case in equation (3.64). The case of zone rings is 


95 

considerably more complicated, because it requires the use of covariances between area- 
mean values as discussed by Colombo (1981a, section 4). For the small size 
compartments (e.g., Ay = 2') that will be considered in the next chapter, the difference 

between the case of rings formed by point measurements and zone rings will be 
neglected. 

Based on the above formulation, it is now possible to examine the particular form 
that the covariance matrices in (3.85) take, in case ring-average values are used for the 
estimation of ATab- The following notation is introduced: 

yj: semi-aperture of terrestrial cap 

y£: semi-aperture of cap at altitude 

Ay T : spacing between rings in the terrestrial caps 

Ay 5 : spacing between rings in the caps at altitude 

yd = yAB : angular separation between stations A and B 
R(- 6371 km): mean-Earth radius 
h: average altitude of satellite observations 

R s = R + h: radius of average sphere where satellite observations refer. 

According to the above, the number of rings in each terrestrial cap is: 

n7=^ + i 

Ay T 

and similarly in each cap at altitude: 

N? = — + 1 

Ay 5 

so that the total number of rings at each site (A or B) is: 

N r = Nj + N 5 


(3.96) 

(3.97) 

(3.98) 


Accordingly, the square symmetric matrix: 
C = C zz (A,A)-Czz(A, B) 
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(3.99) 

appearing in (3.85b), can be arranged as shown in Figure 12, after the observations on 
the Earth and at altitude are collected in two separate groups. 


N 8 — 



Figure 12. Structure of Covariance Matrix C = C ZZ (A, A) - Czz(A, B), and Vector U of 
Equation (3.85). 

Abbreviating by P n (i) the Legendre polynomial P„(cos\|q), where \jq is the spherical 
distance from a cap center to the i-th ring around it, one can easily see that due to the 
equation (3.95c), the elements of C are given by: 

\ 

oo 

C„(i, j)=X WMIl - PJy/JH (a) 

n=2 

oo 

Cl2(iJ)= X <J n(^'| n+2p n( i ) P n(jIl ' P n(Vd)] 0 3 ) j (3.100) 

oo 2 

Cwft j)-2 P„(MI1 - PjVaS «=) 


I 




while, due to (3.95b), the elements of vector U are given by: 


97 


u l(i) - X n + jd n P n (i)[l - P n(Vd)] 

n=2 

Obviously in equations (3. 100) and (3. 101), Pn(Vd). is used to abbreviate P n (cos\j/d). In 
addition \|q is given by: 

(i - 1)A\|/ T (terrestrial data) 

(i-ljAxj/ 5 (data at altitude) . (3.102) 

Finally, the a-priori error variance of the difference ATab (he., the error variance before 
the introduction of the measurements), is twice the quantity: 



(a) 


(b) 


(3.101) 


Ctt(A, A) - Ctt(A, B) = X - P n(Vd)] , 

n=2 ' n+1 ' (3.103) 


as it can be easily deduced from (3.95a). 


The last matrix which needs to be defined so that equation (3.84) can be evaluated, 
is the noise covariance matrix D(A, A) appearing in (3.85b). Assuming no correlation 
between the errors of the terrestrial data and the data at altitude, matrix D(A, A) will be 
block diagonal with one block having dimension N? while the other N, . If On and o£ are 
the degree variances of the error of a point gravity disturbance measurement on the 
Earth’s surface and at altitude respectively, then the element d(i,j) of D(A, A) 
representing the error correlation between the i-th and j-th ring-average value is given by: 


d(i, j) = X p n(i)Pn(j)< 

n=2 


oI\ 

o$) 


(terrestrial data) 
(data at altitude) 


(3.104) 


98 

The above formulation for lsc using ring-averages treats the general case where 
observations on the surface of the Earth and at altitude are used simultaneously for the 
estimation of ATab- The problem of estimating ATab from either group of data alone, 
can easily be treated by forming only the pertinent part of the matrices and vectors 
previously defined. 

3.4 Covariance Models for Signal and Noise 

The implementation of the formulation given in sections 3.1 and 3.3, for the 
estimation of the global mean square error of the disturbing potential differences, requires 
appropriate covariance models to be established that describe, in a global average sense, the 
properties of the gravitational signal and the noise contained in the measurements. The 
various covariance models considered for this purpose are discussed next. 

3.4.1 Global Covariance Models for the Gravity Anomaly 

It is well known (Moritz, 1972) that a global covariance model for the gravity 
anomaly carries equivalent information with a corresponding model for the disturbing 
potential. Hence, covariance models for the gravity anomaly will be considered next, and 
the analytical relations will be given to derive the corresponding models for the disturbing 
potential or other linear functionals of it (e.g. the gravity disturbance). The spherical 
harmonic expansion of the disturbing potential is written as: 


ffl” X (wMe, >.) 

n=2 m=-n (3.105) 

where the notation is consistent with that used in equation (2.52). In an equivalent form 
(3. 105) may be written as: 

ifr, 0, X) = ^ £ (J) n+1 X C^Y 4e, X) . 

n=2 n.=-n (3.106) 


Consistent with the form (3.106), the spherical harmonic expansion of the gravity 
anomaly is: 
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Ag(r,9,X) = 2MX( n ->)(f) n * 2 Z fVMeA) 

n _2 m=-n 

and of the gravity disturbance: 

oo ^ 

6g(r,e.>.) = SMX( n + 1 ^)" +2 S C„mY Je, x) . 

a 1 n=2 m=-n 


(3.107) 


(3.108) 


The following values will be adopted here for the scaling factor "a" and the geocentric 
gravitational constant GM in the above and in the following equations: 

a = 6378137. m \ 

GM = 3986004.36 x 10 8 m 3 /s 2 l (3.109) 


The surface spherical harmonic of the disturbing potential, referring to the surface of the 
sphere of radius: 

r = R = 6371000. m (3.110) 

is given by: 


T„(r = R, 0, A.) = T n = ^j^(if') n+1 X CnmYnmlO, A.) 


(3.111) 


while from (3.107) and (3.108) it can be easily seen that the corresponding surface 
harmonics of Ag and 5g are: 


Ag n (r = R, 0, A) = Ag n = H^-T, 
8g n (r = R, 0,A) = 8g n =n±l-T n 


(3.112) 

(3.113) 


The degree variances (Heiskanen and Moritz, 1967, p. 259) of T, Ag and 5g > referring to 
the surface of the sphere of radius R, are given by: 
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(a) ' 


C n = M(Agn) = (^) 2 kn 

1 

(b) 

(3.114) 

d n = M(6gn) = ( IL ± J -Fk n 

(0 j 


The spatial covariance function of the disturbing potential is given by (Moritz, 1972, 
equation 7-26): 

K(P, Q)= Xk n (i^) n+1 Pn(cos Ypq) 


(3.115) 

while, due to the law of propagation of covariances (ibid, p. 97), the corresponding 
covariance functions of Ag and 5g are: 

oo ^ 

C(P, Q) = X C n(^f + P n(cOS VPQ) 


(3.116) 

and: 



D(P,Q) = £ d^pVjcos Vpq ) . 


(3.117) 

In the notation of section 3.3 one has: 



c t .t(p.Q)=k(p.Q) 


(3.118) 

c 5 g.8 S (P.Q) sD <P.Q) 


(3.119) 


while by the law of propagation of covariance: 
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Ct.5«(P. Q) = - ^-[Ct.t(P. Q)] 

OTQ 




(3.120) 


Based on equations (3.115) through (3.120), all covariances needed in section 3.3 
can be defined. Due to the relations (3.1 14), a model for c„ is necessary and sufficient to 
enable the evaluation of all the above covariances. Up to a certain maximum degree N max 
(which at present equals 360), the anomaly degree variances c„ may be obtained from the 
harmonic coefficients of the disturbing potential, estimated by combining a lower-degree 
set of satellite-derived harmonics with the harmonics implied by a global set of gravity 
anomalies. A state-of-the-art global geopotential solution of this kind is the OSU89B 
model (Rapp and Pavlis, 1990), which is complete to degree and order 360. Hence, up to 
Nmax = 360, c n may be defined by: 



2 £ n £ 360 


(3.121) 


where Qm, are the fully-normalized unitless spherical harmonic coefficients of the 
OSU89B model (even zonals are remainders after removing the ellipsoidal reference 
values). Since the disturbing potential is not band-limited to N max = 360, one needs to 
adopt a model that provides the anomaly degree variances above this degree. The 
analytical form of degree variance models developed for this purpose, is usually selected 
in such a way that enables closed expressions to be derived that provide the covariance 
between various gravimetric quantities, facilitating in this manner the practical 
implementation of these models (Tscheming and Rapp, 1974). The parameters of these 
models are then estimated in a least-squares adjustment where "observed" degree variances 
are in essence the values given by (3.121) from the analysis of satellite and terrestrial data. 
In the following discussion two such models are considered, both of the same analytical 
form (ibid, 1974): 


A(n - 1) 

C " (n-2Xn + B) 


s n +2 


n > 2 


(3.122) 


where A, B and s are parameters of the model; s is the ratio: 
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(3.123) 


where R B is the radius of the embedded (Bjerhammar) sphere. The first set of parameters 
(A, B, s) are those estimated by Tscheming and Rapp (ibid, p. 22) so that the model 
(3.122) is given by: 

0 ^ = 7.5 mgal 2 

cO) = ip^i^il(o.999617) n+2 mgal 2 ; n>2 (3.124) 

The second set of parameters (A, B, s) was estimated by Jekeli (1978) and imply the 
model: 


c^ = 7.5 mgal 2 

<£) = 3 ^ 3 - 3 ^ 8 | n 2 ^l (0.9988961) n+2 mgal 2 ; n > 2 

In Figure 13, the anomaly degree variances implied by the OSU89B model (to 
Nmax = 360) and by the models C$0 = 1, 2) above, are shown. For comparison 
purposes, the anomaly degree variances implied by Kaula’s rule of thumb are also shown. 
These were computed by (Rapp and Pavlis, 1990): 

mgal 2 . (3126) 

Note that expressions (3.121), (3.124) - (3.126), all provide the degree variances 
referring to the surface of the sphere of radius R = 6371 km. As it can be seen from 
Figure 13, from the three models c$, c ( $, c ( „ ) the model c (2) implies the fastest decay of 
the gravity anomaly spectrum at the higher degrees. In addition, one observes a very good 
agreement between the model c (2) and the ("observed") degree variances implied by 
OSU89B in the range 180 < n < 360. Although the spectral properties of the disturbing^ 


(3.125) 


LOG 1 0 (DEGREE) 


Figure 13. Anomaly Degree Variances Implied by OSU89B and the Models c™, c ( n , cv 
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potential above degree 360 are to a large extent a matter of assumption, due to the lack of 
detailed data on a global basis, the overall performance of a degree variance model may be 
judged by comparing the variance of certain gravimetric quantities, as implied by the 
model, to corresponding values obtained from the analysis of available measurements. 
For this purpose, the variance of the (point) gravity anomaly Co, and the variance of the 
horizontal gradient of the gravity anomaly Goh are often used (Moritz, 1980, section 22). 
Jekeli (1978) has shown that: 


Goh * ~ 

L n=2 


where the degree variance of the vertical gradient, g n is: 
(3Ag\ _ (n + 2 f r 


gn = 


l ar /n R 2 


(3.127) 


(3.128) 


To calculate Co and Goh, four models for c n are considered as follows: 


Cn(A) S cS;> 


III 

sr 

o 


c«(C) = j 

fc,i(OSU89B) 

14" 

2 < n < 360 
360 < n < oo 

C n(D) = | 

k(OSU89B) 

l<4 2) 

2 £ n £ 360 
360 < n < <*> \ 


(3.129) 


Using these models and approximating infinity by n = 36000, Co and Goh have been 
evaluated. The values obtained are given in Table 4. 
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Table 4. Gravity Anomaly (Co) and Horizontal Anomaly Gradient (Goh) Variances, 
Implied by Different Anomaly Degree Variance Models. 


Model 

Co(mgal 2 ) 

Goh(E.U. 2 ) 

A 

1795.0 

3542.2 

B 

1106.0 

339.3 

C 

1437.9 

3533.2 

D 

1035.8 

339.3 


Based on extensive analysis of gravity anomalies implied by satellite altimetry, as well as 
of terrestrial measurements, Rapp (1985) suggested 1100 mgal 2 as an estimate for Co. 
The value of Goh is much more difficult to estimate, since the horizontal gradient is a 
signal of high-frequency content and the detailed data necessary for its determination are 
not currently available on a global basis. A value of 300 EU 2 (1 Eotvos Unit (EU) = 
0.1 mgal/km = 10* 9 s 2 ) for Goh is a compromise between various estimates obtained from 
regional data analyses (Robbins, 1985). Based on these considerations, the composite 
model c n (D) defined in (3.129), will be used in the numerical analysis to be presented in 
the next chapter. 

3.4.2 Covariance Models for the Measurement Noise 

The estimation of realistic error properties for the data poses an even more difficult 
problem than the estimation of the global properties of the gravitational signal contained in 
them. Provided multiple independent samples of die measurements are available, one may 
estimate an empirical error covariance function for the data by analyzing the differences 
between different samples, as it was done by Weber and Wenzel (1982) for the gravity 
anomaly data in Europe. In the absence of multiple independent measurements for an 
observable, one usually employs an analytical error covariance model to describe its error 
properties. The variance and the correlation length of such a model are assigned values 
which are representative of the measurement process under consideration. In the 
following, three such error covariance models are considered. 



Model 1 : Dirac Impulse 
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The error covariance function here is given by: 

°^ l \v) = cD(v) ; c>° (3.130) 

where the Dirac impulse (or delta function), D(y), is defined on the unit sphere by Jekeli 
(1981): 


4n 


jj d(ypq) f(e P , x P )do=f(eQ,XQ) 


(3.131) 


For the isotropic error covariance models that will be considered here, the n* -degree 
variance c n is given in general by: 


o n = + 1 I o(v)Pn(cosy)sin\|/dv(r 

(3.132) 


This equation yields for the model ori)(\|/), due to (3.130) and (3.131): 

</ n l) = (2n + l)c ; n£0 (3.133) 

which implies that the rms amplitude by degree of the error is constant, equal to c 1 ^ 2 , for 
any degree n. This property is the spherical counterpart of the corresponding spectral 
property of the delta function, 8(x), defined on the real line -°° < x < °° whose Fourier 
transform is one (Papoulis, 1962, p. 36). 

The model aO)(Y), as it can be seen from (3.130) and (3.133), implies infinite 
variance and zero correlation length. In practical applications, where one is always 
considering a finite number of observations (point values or area-means), the assumption 
that observation errors are uncorrelated and that all observations have the same error 
variance, is often employed. However, as seen from the above, such assumptions, when 
extended to the case of infinite number of point measurements covering the entire sphere, 
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lead to unrealistic spectral properties of the noise. For this reason alternative error 
covariance models are considered next, which do not suffer from such disadvantages. 

Model 2 : Gauss-Markov First-Order 

The error covariance function here is given by (Rummel, 1980, p. 44): 

d 2 X V) = ce'* , v ; c>0, X>0 . (3.134) 

The parameters c, X are related to the variance (mo) and the correlation length (xj/ 0 ) by: 
c = m$ (a) \ 

J (3.135) 

* = -L/n2 (b) 

V s / 


The degree variances corresponding to 0( 2 Ky) can be obtained from the following 
recurrence relation: 


q(2) _ 2n + 1 
n 2n - 3 


^ 2 + (n -If Jz) 
X 2 +(n + l J 2 n ' 2 


n£2 


with starting values: 

_£ 1 + e - ^ 1 

2 X 2 +l 




1 - e-** 
X 2 + 4 



(3.136) 


(3.137) 


The detailed derivation of this recurrence is given in Appendix C. In case one postulates 
the absence of zeroeth- and first-degree harmonics from o^0|0> the resulting covariance 
model is given by: 
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Given the correlation length xj/ 0 , the parameter X for this model can be determined by 
solving iteratively the equation: 

1 - 2e-' u *' + -L I +c~ x * + i-^^(2cosv c - 1) = 0 

2 X 2 +l 2 H 2 + 4 l (3.139) 


with initial value for X , the one given in (3.135b). The parameter c can then be 
determined, given the variance m?. by: 


c = rr£ 


i _J_ 1+e-** . 3. 1 - 

1 ^ n 


X 2 + 1 


X 2 + 4J 


(3.140) 


Obviously, the degree variances o(, 2 ) are identical to o^for n £ 2, while 

a?)-cp>-a 


Model 3 : Reciprocal Distance 

The error covariance function here is given by: 

o(3^y) = c{i - X)(l - 2Xcos\|r + X 2 ) ; c>0,0<X<l . (3.141) 

Considering the notation of Figure 3, it can be easily seen that o< 3 >(y) can be written in the 
form: 

0 (J) (V ) „ , (3.142) 

where: 

r = ^ (a) 

X 

l = r(l - 2Xcosy + X 2 )^ 2 (b) 



(3.143) 
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and the form (3.142) explains the characterization "reciprocal distance". Given the 
variance mo and the correlation length xj/ 0 , it is easy to see that the parameter c is given by: 

c = mg , (3.144) 

while X is obtained after solving iteratively the equation: 

(1 - 2Xcos\|/c + X 2 ) 1 / 2 - 2(1 - X) = 0 (3. 145) 


with initial value for X given by: 


X (o) =l -Vising 


(3.146) 


The degree variances are given by the closed form expression (Sjoberg and Fan, 
1986): 


ofMl -X)X" 


(3.147) 


The assumption of absence of zeroeth- and first-degree harmonics from a^ 3 \\|/), yields the 
model (ibid, 1986): 


</ 3 \y) = c(l-X)[(l - 2Xcosxj/ + X 2 ) ^ - 1 - Xcosy] 


(3.148) 


The parameter X for the model c? 3 \xj/) is obtained by solving iteratively the equation: 


2(1 -x)[(i - 2 X 008 x 1 /° + X 2 ) - 1 - Xcosx|/°] - X 2 = 0 , 


(3.149) 


given Xj/ 0 , and using as initial value for X the one given in (3.146). Once X has been 
determined, c is given by: 


c = mo X' 2 


(3.150) 



Again, for n £ 2, while o(, 3 ^ = = 0. 
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In order to demonstrate the different characteristics of the models 
o (2 W) and o^ 3 \\}/), the parameters c and X of these models were evaluated for a given 
value of the variance m^ (taken here to be m© = 1 mgal 2 ) and different choices of the 
correlation length (taken as 0°1, 0°2 and 0°5). The results are given in Table 5. For 

the model c^ 3 \\|/) the radial distance r to the external point P (see Figure 3) from which the 
distance l is reckoned is also given. The radius R in this case (see equation 3.142) was 
taken to be 6371 km. 

Table 5. Parameters Associated with the Error Covariance Models c/ 2 \\|/) an d </ 3 \y), for 
mo = 1 mgal 2 and \\f° = 0^1, 0°2 and 0°5. 


Model 


0?1 

0?2 

0°5 

o (2 W) 

H 

0.397 1 368 1 487 3D+03 
1.000012680788 

0.198557 5090 1 3D+03 
1.000050727516 

0.793924826289D+02 

1.000317237947 

ct ( 3 W) 

X 

c 

r(m) 

0.998990129118 

1.002022805406 

6377440.4 

0.997975830838 

1.004060663364 

6383922.1 

0.994906095858 

1.010266183954 

6403619.4 


As expected, for a given value of the variance mo, increase of the correlation length 
causes decrease of the parameter X for both models c/ 2 \y) and c? 3 \y) (see equations 
3.135b and 3.146). As \|^ increases, more power of the error spectrum is shifted to the 
lower degree harmonics, as it is illustrated in Figures 14a and 14b where the degree 
variances and are plotted for \\r c = 0°1 and \\f c = 0°2 respectively. For the 

reciprocal distance model c? 3 \\y), increase of vj/ 0 corresponds also to increase of the radial 
distance r, as it can be seen from Table 5. 




















CHAPTER IV 


NUMERICAL ANALYSIS 

Based on the formulation given in Chapters II and III, a number of numerical 
experiments were performed in order to assess the accuracy of geopotential differences 
estimated from the gravimetric information contained in a global spherical harmonic 
expansion and in the gravity disturbance data in the terrestrial and at altitude caps. The 
results of these computations are presented in order of increased complexity of the 
estimation scheme employed, starting with the geopotential differences estimated on the 
basis of a global geopotential model alone and concluding with the addition of gravity 
disturbance data on the surface of the Earth as well as at altitude. 

4.1 Disturbing Potential Difference Estimated from Current and Future 

Global Geopotential Solutions Only 

Three global geopotential solutions are considered here. These are designated: 

(a) OSU89B 

(b) TOPEX 

(c) GPB 

(a) OSU89B (Rapp and Pavlis, 1990) represents a state-of-the-art high-degree global 
model, available in the absence of the missions discussed in section 1.3. The estimated 
anomaly error degree variances associated with this model are shown in Figure 13. A 
number of comparisons with independent data reported by Rapp and Pavlis (ibid), indicate 
that the error estimates of the coefficients do provide a realistic assessment of the quality of 
the model. As it can be seen from Figure 13 the spectrum of the anomaly error starts 
exceeding the spectrum of the anomaly signal around degree 250. As the degree of the 
coefficients increases, validation of their error estimates becomes more difficult because 
precise independent information is not available on a global basis with high enough 
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resolution to enable comparisons. In ocean areas comparisons with GEOSAT-implied 
undulations indicated that, although the signal sinks below the noise above degree 250, the 
coefficients above this degree do contain meaningful information (ibid, p. 21,903, Table 
9). In this study for the implementation of least-squares collocation, the OSU89B 
coefficients above degree Nmax = 250 will not be admitted as part of a reference model, 
with the understanding that the pessimistic error estimates above this degree may be due to 
shortcomings in the error assessment of the model, induced by factors such as the neglect 
of correlations between the 30' x 30' mean anomalies used in its development. Such 
neglect is at present necessitated by the inability to computationally handle a more 
appropriate error modeling for the 259200 mean anomalies involved in the model's 
development. 

(b) TQPEX represents the geopotential model that is anticipated to become available from 
the analysis of the GPS tracking data alone on the TOPEX/Poseidon altimeter satellite 
(Pavlis, E., et. al., 1990). Such a model will extend to degree 50 and its expected errors 
were estimated in a simulation study where the GPS tracking data were used in a dynamic 
solution for the estimation of the geopotential spectrum, satellite orbital parameters as well 
as parameters related to air-drag etc. (ibid, 1990). 

(c) GPB represents the geopotential model anticipated to become available from the 
analysis of the GPS tracking data on Gravity Probe-B spacecraft The model from such a 
mission will extend to degree 55, and the expected errors of its coefficients were obtained 
in the simulation study by Pavlis, E., et al. (ibid) in a similar manner as for TOPEX. 

The error anomaly degree variances for OSU89B (2 < n < 60), TOPEX and GPB 
are shown in Figure 15, along with the (signal) anomaly degree variances implied by 
OSU89B (2 £ n <> 60). As it can be seen from this figure, TOPEX is expected to improve 
the current knowledge of the spectrum at degree 10 by about 3 orders of magnitude while 
the improvement for GPB at the same degree is by about 4 orders of magnitude. 
However, since TOPEX and GPB models lack the detailed terrestrial (and altimetry- 
implied) gravity information included in OSU89B, their error spectra yield poorer values 
than the error spectrum of OSU89B, after degrees ~ 25 and - 45 respectively. 
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Figure 15. Anomaly Degree Variances for OSU89B and Error Anomaly Degree Variances 
for OSU89B, TOPEX and GPB. 

Since the information based on which TOPEX and GPB are developed is 
independent of the information used to develop OSU89B, one can safely assume that if 
OSU89B is combined (in a least-squares sense) with either of the two models, the 
resulting solution will perform at least as good as the best of the contributing solutions 
performs at any given degree. In that sense, the error spectrum of a combined solution of 
the type TOPEX/OSU89B may be approximated by the error spectrum of TOPEX up to 
degree 25 and the error spectrum of OSU89B from 26 to 360. Similarly, a combined 
GPB/OSU89B model's error spectrum can be approximated by the error spectrum of GPB 
up to degree 45 and the error spectrum of OSU89B from 46 to 360. Based on these 
considerations, the total (commission plus omission) error in geopotential (Yd = 0°) and 
geopotential difference, implied by the above models, truncated at various degrees was 
computed using equations (3.53), (3.54) and (3.55). The results are given in Table 6. In 
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all numerical results to be presented hereon the radius of the reference sphere to which the 
estimated errors refer is R = 6371 km, and the model c n (D) of section 3.4.1 is used to 
represent the "true" spectrum of the gravity anomaly. 

Table 6. Geopotential and Geopotential Difference Errors from Current and Future 
Global Spherical Harmonic Models. 



eAT (kgal cm) 

Model 

Nmax 

Vd = 0° 

T 


1(F 

3(F 

180° 

OSU89B 

25 

207.7 



291.6 

294.0 

290.7 

(89B) 

45 

154.6 

Iffli 

1 

221.7 

217.9 

217.1 


250 

60.6 

77.0 

1 

86.2 

85.6 

85.9 


25 

206.8 


310.5 

290.2 

292.7 

289.4 

II 

50 

161.3 

BlilB 

236.4 

228.0 

228.5 

228.9 

GPB 

25 



HRQ| 

289.8 

292.5 

EH 

(G) 

50 

142.5 

HM 

Wnmm 

205.6 

202.1 

EE99 

T/89B 

25/250 

57.4 

76.5 

83.0 

81.2 

81.1 

81.2 

G/89B 

45/250 

52.2 

74.2 

74.7 

73.9 

73.8 

73.8 


As it can be seen from Table 6, the current knowledge of the geopotential spectrum 
implies absolute geopotential values (yd = 0°) accurate to about 61 kgal cm, while for 
station separation = 10° (~1 1 10 km) the error in AT is about 86 kgal cm. The future 
models alone, (i.e. N max = 50), despite the great improvements they promise for the very 
low-degree harmonics, fail to achieve the accuracy obtainable by the currently available 
high-degree expansion OSU89B, due to the large omission error above the maximum 
degree that these models can resolve. However, when the high accuracy long-wavelength 
information provided by GPB, is augmented by the shorter-wavelength information 
contained in the high-degree model, an improvement by about 14% over the attainable 
accuracies using OSU89B alone (up to N m ax = 250), is achieved. The corresponding 
improvement for the combination TOPEX/OSU89B is substantially lower (about 5%) as it 
is expected (see Figure 15). Apart of the higher accuracies achievable by the combination 
GPB/OSU89B, an additional factor that makes this combination preferable over the 
TOPEX/OSU89B, is that the higher resolution of GPB enables better control over the 
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systematic errors in the combined solution, caused by the errors in the gravity anomalies 
(used in the development of OSU89B), arising from vertical datum inconsistencies. As 
mentioned in section 1.1, Laskowski's (1983) study indicated that such errors are likely to 
contaminate the geopotential's spectrum up to maximum degree about 60. 

4.2 Introduction of Terrestrial Gravity Disturbance Measurements 

The estimation scheme of the previous section, whereby the geopotential 
differences were obtained on the basis of the global information contained in a geopotential 
model only, is now augmented by introducing detailed local gravimetric information in the 
form of gravity disturbances inside the caps surrounding the computation points. The 
characteristics and sensitivities of this estimation scheme are considerably more 
complicated to assess, than of the scheme of section 4.1, because here a number of inter- 
related factors affect the quality of the final result. Hence, to examine the individual 
contribution and importance of different aspects of the observation/estimation setup to the 
quality of the results, the following discussion is divided into three parts, addressing the 
choice of cap integration kernel, the cap size and data density, and the data accuracy 
respectively. In addition, since the error estimation of the observational setup considered 
here, can be performed using either truncation theory principles or least-squares 
collocation, a comparison of the error estimates obtained by each technique is also given in 
section 4.2.4. 

4.2.1 The Choice of Cap Integration Kernel 

The two integration kernels Hi(\jO (unmodified Hotine) and H 2 (y) (Meissl's 
modification) are examined here. The propagated errors in geopotential (yd = 0°) and 
geopotential difference implied by the use of these kernels are given on Table 7, for two 
cap sizes (\|/q = 2°, 4°) and two choices of reference geopotential model (OSU89B and the 

combination GPB to degree 45 augmented by OSU89B from 46 to 250). 
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Table 7. Geopotential and Geopotential Difference Errors Implied by Different Choices of 
Cap Integration Kernel. (All Errors are in kgal cm). 




OSU89B 

GPB/OSU89B 

ma 

Vd 

Hi(v) 

h 2 (v) 

Hi(v) 

h 2 (v) 


OP 

19.4 

19.8 

16.7 

msm 

2 ? 

5 > 

24.4 

24.3 

22.6 

EH 


10P 

27.2 

28.5 

23.4 

EH 


30P 

27.5 

28.0 

23.6 

12.3 


OP 

18.6 

mm 

16.3 

MEM 

4? 

5 * 

23.5 

HI 

19.9 

EH 


10P 

24.8 

EH 

21.8 

HI 


30° 

26.3 

19.8 | 

23.1 

■EH 


T 

Vo : terrestrial cap size 

Yd • angular separation between stations 

AyT = 4' 

Gravity Disturbance Error : a^ 3 \y) ; mo = 2 mgal, tj/ 0 = 0® 1 
Nmax = 250 


In the evaluation of the errors given in Table 7 (as well as in all subsequent 
applications of truncation theory) the sampling error spectrum s„, which was introduced 
but not specified in equation (3.51a), was evaluated from the quartic expression developed 
by C. Jekeli based on numerical experiments performed by Colombo (1981a, equation 
3.10), i.e. 


Sn_f 1'N') dn 

where: 

= 16-1957^) + 30.34506](-JL) + 40.29588^ 


(4.1) 


(4.2) 


and N is the Nyquist frequency implied by the sampling interval (equation 3.21), while d n 
is the signal degree variance of the gravity disturbance. 


Examining now the results shown in Table 7, the following comments can be made: 
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1. By comparing the results for the case of OSU89B and the two different cap sizes it can 
be seen that the kernel modification yields improved accuracy estimates by about 25% for 
the larger cap size (4°), but degrades slightly the results for the smaller cap size. In 
contrast, in the case of the combined model GPB/OSU89B, the kernel modification yields 
a substantial (~ 48%) reduction of the error for the smaller cap size and a less pronounced 
(~ 39%) for the larger one. This is due to the fact that Meissl's modification, in its attempt 
to accelerate the rate of convergence of the truncation coefficients, shifts part of the power 
of the higher-degree coefficients to the lower-degree ones (see also (Jekeli, 1980). 
Accordingly although the omission error decreases with the introduction of Meissl's 
modification, the commission error increases, and that explains the slight degradation of 
the results when OSU89B is used as reference model. In case the combined model is 
used, the lower-degree harmonics of this model are so accurately determined from GPB, 
that the increase of the commission error is overwhelmed by the concurrent decrease of the 
omission one, and that accounts for the great improvement achieved by the use of H 2 O 1 O 
instead of Hi(\j/) here. 

The increase of the cap size to 4° implies that additional information concerning the 
long-wavelengths is contributed by the cap measurements, as compared to the case 
\|fj = 2°. This counterbalances the increased commission error in the case of OSU89B 
(going from Hi(y) to H 2 O 1 O), but, since the cap measurements are not error-free, reduces 
the 48% improvement in the case of the combined model, to 39% (increased cap size 
implies increased propagated error). 

2. A comparison of the results for the case (Vo = 2°, H 2 O 1 O, GPB/OSU89B) in Table 7, 
with the results given in Table 6 for GPB/OSU89B (Nmax = 250), shows the substantial 
improvement achieved by the incorporation of terrestrial gravity disturbances around the 
computation points. Such comparison implies that the present accuracies of the higher- 
degree harmonics, as well as the resolution that current high-degree expansions achieve, 
are only capable of geopotential difference determination at the accuracy level of about 75 
kgal cm. If one strives for a 10 kgal cm accuracy level, then, not only terrestrial 
gravimetry, but the high-accuracy long-wavelength geopotential determination anticipated 
from the GP-B mission are needed. As it can be seen from Table 7, the present 
knowledge of the long-wavelength (2 to 45) part of the geopotential spectrum, can only 
support geopotential difference determination at the 20 to 30 kgal cm level for the cap sizes 
that were considered. 


4.2.2 The Effect of Cap Size and Data Density 
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The dependence of the accuracy of the estimated geopotential differences on the 
extent (cap size) and density of the gravity disturbance measurements in the caps 
surrounding the computation points is very important, since precise gravimetric surveys 
required to gather the data are rather costly undertakings. In order to examine these 
aspects, the error of the estimated geopotential (\j/{j = 0°) and geopotential difference was 
evaluated for three cap sizes (\|/$ = 1°, 2° and 4°) and three grid spacings (A\|/ T = 2', 6' and 
10'), with the combined GPB/OSU89B model to Nmax = 250 used in all cases as the 
reference and with H 2 (y) as the cap integration kernel. In all cases a reciprocal distance 
error covariance model with variance 4 mgal 2 and correlation length 0°1 was used. The 
results are given in Table 8. 

Table 8. Influence of Data Extent (Cap Size) and Data Density on the Error of Estimated 
Geopotential and Geopotential Difference. (All Errors are in kgal cm). 
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As it can be seen from Table 8 the grid spacing (Ay 7 ) has only a minor effect on the 
quality of the result. As one passes from the very dense sampling (Ay T = 2'), to the 
coarser one (A\j/ T = 10’), a degradation of the results by only about 1 kgal cm is observed. 
This is due to the fact that the sampling error, for sampling intervals smaller than ~ 10’, 
has the smallest contribution to the total error of the geopotential and the geopotential 
difference. For the cap sizes under consideration (up to ~ 5°) and the reference high- 
degree spherical harmonic models used here, the two major sources of error contributing 
to the total error budget are the propagated error of the cap data and the commission error 
of the model. The little dependence of the total error on the sampling interval (for Ay T < 
10’), is rather fortunate because it implies cost savings in the data acquisition process. 

As far as the cap size is concerned, as it can be seen from Table 8, Vo = 2° yields 
the best results for the particular selection of reference gravity model (and its higher 
degree) and error properties of the data. It is important to notice that due to the 
accumulation of propagated error of the data, the increase of the cap size beyond a certain 
value does not improve but degrades the results. Although this appears strange at first 
glance (introduction of additional information should not worsen the determination of T 
and AT), it is well explained if one takes into account that the kernel H 2 (y), used to 
evaluate the errors in Table 8, is designed based on deterministic principles with no 
explicit requirement to minimize the total error of the geopotential. This aspect will be 
reconsidered in section 4.2.4, since it constitutes a major difference between integral 
formulae techniques using deterministic kernels and least-squares collocation. 

4.2.3 The Effect of the Error Properties of the Gravity Disturbance Data 

In view of the difficulties encountered in estimating realistic error properties for the 
gravity disturbance data inside the caps, it is important to examine the effect of various 
assumptions concerning the behavior of the data errors, on the errors of the estimated 
quantities (T and AT). For this purpose, the two error covariance models 
o^y) and o 3 (y) were used, with different variances (mo =1,4 mgals 2 ) and correlation 
lengths (Vo = 0°1, 0°5) for the estimation of eTtyd = 0°) and eAT. The results are given 
in Table 9. In all cases, the cap size was 2°, the sampling interval was 4’, the geopotential 
model GPB/OSU89B to N max = 250 was used as reference model and H^y) was used as 
integration kernel. Examining the results in Table 9, the following comments can be 
made: 
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Table 9. Influence of the Error Properties of the Gravity Disturbance Data on the Error of 
Estimated Geopotential and Geopotential Difference. (All Errors are in kgal cm). 
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1. In all cases considered the Gauss-Markov model o^ 1 2 \x|f) yields more optimistic results 
than the reciprocal distance model o^ 3 \y). Observing Figure 14a, one can see that the 
first-order Gauss-Markov model contains more power of the error below harmonic degree 
~ 150, than the reciprocal distance error model does, for the same degree range. Since the 
cap data error contribution at the lower degrees is greatly attenuated, due to the very small 
magnitude of the coefficients XQkn (see equation (3.50b)) at these degrees, for a given 
variance m§, the Gauss-Markov model o^ 2 \x|/) implies a smaller propagated error, than 
the reciprocal distance does. In addition, the difference of the errors implied by the two 
models becomes smaller as the correlation length increases, because the increase of y c 
shifts more power of the error ci^ at lower degrees as it can be seen from Figure 14b. 
For example, for m3 = 1 mgal 2 and xj/ 0 = 0°1, a^ 3 \\ji) yields an error ~ 8% larger than 
a^ 2 \xy), while the difference is only ~ 3% for xj/ 0 = 0®5 (for the same error variance). 

2. The error of T and AT is affected more by changes of the error correlation length, than 

by changes of the error variance. For example, for = 4 mgal 2 , the error nearly doubles 
as the correlation length increases from 0°1 to 0°5, while for x|/ c = 0°1, increasing the 
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variance from 1 to 4 mgal 2 causes only an increase of the error of approximately 26%. 
Note also that for the larger variance, an increase of the correlation length causes a quite 
larger increase of eT and eAT, than for the smaller variance. 

According to the above, larger but less correlated errors in the cap measurements 
can be tolerated more than smaller but heavily correlated ones (e.g., systematic errors in 
the cap data), as far as the determination of the geopotential differences is concerned. 

4.2.4 Comparison Between Error Estimates from Truncation Theory and 
from Least-Squares Collocation 

The two different techniques for assessing the global rms error of the geopotential 
differences, obtained on the basis of a global spherical harmonic model and gravity 
disturbance data in the terrestrial caps, are intercompared in this section. For this purpose, 
the parameters given in Table 10 were used in the implementation of truncation theory, and 
the error eAT was evaluated with both techniques, for two cap sizes (V5 = 2°, 4°), using 
the two different error covariance models, for three station separations (\|/d = 5°, 10°, 30°). 


Table 10. Comparison Between Error Estimates of Geopotential Differences Obtained 
from Truncation Theory and Least-Squares Collocation (All Errors are in 
kgal cm). 
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Examining the results given in Table 10, the following observations can be made: 
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1. The use of least-squares collocation (with ring averages) always yields smaller error 
eAT than the use of the integral formula. This is as expected, since lsc is by definition the 
optimal linear estimator, in the sense that it minimizes the rms error of the estimates, 
among all the other linear estimators, one of which is the integral of Hotine. 

2. Increase of the cap size can never cause increase of eAT in the case of lsc, no matter 
what the error properties of the data are. This is not the case when the integral formula is 
used, as it was also seen in section 4.2.2, due to the accumulation of propagated error 
caused by the increase of the number of measurements. 

3. The reciprocal distance error covariance model yields more pessimistic results than the 
Gauss-Markov model in the case of lsc, as it was also the case when the truncation theory 
was used to estimate the error eAT. 

4. The error estimates from the two techniques are generally in good agreement, the 
largest difference being on the order of 3 kgal cm (~ 21%). 

Finally, summarizing the results of the previous sections, the error analysis 
performed here indicates that with the combination GPB/OSU89B model, used to 
N m ax = 250, and with gravity disturbance data on a 2° cap, accurate to 2 mgals, 
geopotential differences accurate to ~ 12 kgal cm are achievable for stations separated by 
30° (about 3300 km). If the error of the cap data reduces to 1 mgal, the corresponding 
eAT becomes about 8.5 kgal cm, which is a comparable figure with the current accuracies 
of geocentric positions estimated using satellite techniques. 

4.3 Introduction of Gravity Disturbance Data at Altitude 

The GPS/GPB SST configuration is considered here, as the observational system 
that provides the gravity disturbance data at altitude. The long-wavelength reference 
geopotential model used in the least- squares collocation (with ring averages) error analysis 
that is presented next, is assumed to be the product of a global dynamic solution, based on 
smoothed values of the SST data, as was discussed in section 2.4. The maximum degree 
of this model is taken initially equal to 20 for two reasons. First, to avoid the systematic 
errors on the gravity disturbance at altitude, arising from the residual centrifugal 
acceleration term 5Ri 0 as explained in section 2.4 (see Table 3), and second, to permit a 
higher signal-to-noise ratio for the data at altitude, as opposed to higher degrees of 
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truncation. Two different caps sizes (Yo = 5°, 10°) were considered for the data at altitude, 
and for each case two different ring spacings were used (Ay s = 30', 60'). Data density of 
one observation per square degree at altitude, implies an integration interval of about 13 
seconds, for the polar 10-day repeat orbit of GPB, based on an approximate formula given 
by Jekeli and Rapp (1980). Considering also terrestrial gravity disturbance data in caps of 
semi-aperature 2° and ring spacing 6', the error estimates eAT for station separation 
Yd = 30° were computed, using also different assumptions for the error properties of the 
data at altitude (rpg = 0.4, 0.2 mgal). The results are given in Table 11. 

Table 11. Geopotential Difference Errors Implied by the Combination of Gravity 

Disturbance Data on the Earth and at Altitude. (All Errors are in kgal cm). 
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As it can be seen from Table 11, the reciprocal distance error covariance model 
o^ 3 \y) yields again the most pessimistic results. This can be attributed to the fact that the 
downward continuation process involving the data at altitude, amplifies the error power in 
the degree range 150-750, in which range the reciprocal distance error model has more 
power concentrated than the first-order Gauss-Markov model (see Figure 14a). Although 
the downward continuation amplifies the error power throughout the band, the lower- 
degree part is dominated in the above estimation scheme by the reference model, while the 
very high part of the error spectrum has small effect on the estimated geopotential 
differences. 
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Overall, the error estimates given in Table 1 1 are of comparative magnitude with 
those obtained in the absence of the data at altitude. In order to enable a more systematic 
comparison between the two cases, in Table 12 the error estimates eAT obtained on the 
basis of terrestrial data only, and various reference geopotential models, are summarized. 
These were computed using least-squares collocation with ring averages to enable a fair 
comparison with the estimates of Table 1 1 . Comparing now the error eAT from Table 1 1 

for the case N ma x - 45, o^ 3 \\|f), mo = 0.4 mgal, \|/§ = 10°, A\j/ S = 60' (15.7 kgal cm), to 
the corresponding 16.9 kgal cm from Table 12, one can see that an improvement of about 
7% is achievable by the introduction of the data at altitude. Although small, such an 
improvement indicates two things: 

(a) The error properties assigned to the data at altitude according to the reciprocal distance 
error covariance model with variance 0.4 mgal and correlation length 0°1, appear to 
provide a realistic assessment of the quality of them. A larger improvement would imply 
that the error properties assigned to these data are too optimistic and in disagreement with 
the error budget based on which the error estimates of the global geopotential solution 
GPB were derived. 

Table 12. Geopotential Difference Error eAT Implied by Current and Future Reference 
Models and Terrestrial Gravity Disturbance Data in a 2° Cap. 
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(b) The introduction of the data at altitude provide a better configuration for the estimation 
of the geopotential differences, as compared to the use of terrestrial measurements only. 

Using pessimistic error estimates as a guideline, one can conclude that 
measurements from the Gravity Probe B mission (in combination with terrestrial gravity 
disturbance data in 2° caps) are capable of providing geopotential differences over 
intercontinental distances, accurate to about 15 to 17 kgal cm. This estimate is an order of 
magnitude better than what can be achieved using the traditional vertical connections based 
on MSL monitoring (about 150 kgal cm). However, lacking gravitational information in 
the medium frequency band (45 < n < 250), the results from GPB are about 5 kgal cm 
worse that those obtained when GPB is augmented by OSU89B in this spectral band (eAT 
= 12 kgal cm). It should be mentioned though that such augmentation does not guarantee 
that the estimated geopotential differences are totally free of systematic errors arising from 
vertical datum inconsistencies, since the harmonics of OSU89B above degree 45 may be 
contaminated by such errors. To enrich the satellite observations with medium frequency 
gravitational signal, requires either a lower flying altitude, or additional on-board 
instrumentation capable of measuring higher-order gradients of the gravitational potential, 
and in this direction the ARISTOTELES mission may provide a significant contribution. 

Finally, it should be emphasized here that the error estimates evaluated for the 
geopotential differences in this chapter, account for errors of the reference geopotential 
model and the gravity disturbance data only. Additional error contribution to the 
geopotential difference will arise due to random errors of the geocentric positions of the 
stations and due to a systematic offset of the origin of the coordinate system from the 
geocenter, as discussed by Colombo (1980) and Hajela (1983). Random position errors 
at the ± 5 cm level (Smith et al., 1985) have small effect on the estimated geopotential 
differences, while the systematic error due to non-geocentricity of the reference frame can 
have significant effect on the intercontinental connections. 



CHAPTER V 


SUMMARY, CONCLUSIONS AND RECOMMENDATIONS 

The problem of estimation of geopotential differences over intercontinental locations 
was re-examined, in order to assess currently achievable accuracies and future anticipated 
improvements. Accurate estimation of the geopotential differences between points located 
at different continents, imply the unification of the vertical datums established in them, 
which at present are defined based on MSL monitoring and thus are (in general) 
inconsistent due to the presence of the Quasi-stationary Sea Surface Topography. 

A review of the proposed techniques for the unification of vertical datums, in 
conjunction with anticipated future satellite missions and in view of the accuracies 
achievable at present for geocentric positioning, indicated that approaches based on the 
combination of gravitational information with high accuracy geocentric positioning, are 
favored at present and in the near future for practical implementation. In this direction, 
extending the ideas put forward by Colombo (1980), an observational setup was 
proposed, whereby gravity disturbance measurements on the Earth’s surface, in caps 
surrounding the estimation points, are combined with corresponding data in caps directly 
over these points at the altitude of a low orbiting satellite, for the estimation of the 
geopotential difference between the terrestrial stations. The gravity disturbance data at 
altitude are inferred from GPS measurements made from the low orbiter to the high- 
altitude GPS satellites, in a multiple-high-single-low Satellite- to- Satellite Tracking 
configuration. In the absence of actual measurements, the performance of such an 
observation/estimation scheme was evaluated by conducting an error analysis study. 

The mathematical modeling required to relate the primary observables to the 
parameters to be estimated, was studied both for the terrestrial data and the data at altitude. 
Emphasis was placed on the examination of systematic effects and the corresponding 
reductions that need to be applied to the measurements to avoid systematic errors. For the 
gravitational accelerations inferred from SST data, it was discovered that the magnitude of 
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a centrifugal acceleration term (5Ri 0 ) was underestimated by several orders of magnitude 
in the past as a result of an erroneous derivation. The previous formulation implied a 
magnitude of 5Rj 0 about 7 x 10 5 times smaller than the current corrected formulation. It 
was shown in this study that in order to keep the systematic effect arising from 5Rj 0 at the 
20 pgal level, a reference geopotential model complete to degree 20 is required (high-low 
SST configuration). Previous analyses, based on the erroneous formulation, were 
indicating that a reference model complete to degree 4 is adequate to keep the residual 
systematic effect of 5Rj 0 at the 10 p.gal level. For a given noise level (0.4 mgal) of the 
data at altitude, increase of the maximum degree of the reference model, significantly 
affects the ratio of the residual signal to the noise. 

Two different techniques were considered for the estimation of the global mean 
square error of the geopotential differences. Error propagation using truncation theory, as 
applied to Hotine’s integral formula, and the least-squares collocation using ring averages 
as input data. Both techniques are applicable in case observations on the Earth’s surface 
only are involved in the geopotential difference estimation, but only Isc can handle 
efficiently the over-determined case when observations at altitude are added. Alternative 
formulations related to the sampling (or discretion) and the propagated errors arising in the 
truncation theory considerations were derived. These are characterized by the same 
computational requirements as the previous formulation by Christodoulidis (1976), while 
they provide a more consistent interpretation of the underlying physical principles that give 
rise to these errors. In an attempt to apply truncation theory principles for the assessment 
of the contribution of gravitational acceleration data at altitude, to the estimation of 
geopotential differences on the Earth's surface, recurrence relations for the altitude 
generalized truncation coefficients implied by Hotine's kernel were developed for the first 
time. 


Both techniques (truncation theory and lsc) require for their implementation a-priori 
knowledge of the global properties of the signals and the noise involved in the estimation 
and to this end different covariance models were considered and their spectral 
characteristics were compared. In addition, an efficient recurrence relation for the degree 
variances implied by a first-order Gauss-Markov covariance model was developed for the 
first time in this study. 
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For the numerical analysis, three global geopotential solutions were considered as 
reference models. The currently available OSU89B high-degree harmonic expansion, and 
the global models anticipated to become available from GPS tracking data of the 
TOPEX/Poseidon and the Gravity Probe B spacecrafts respectively. Augmentations of the 
latter two models with higher-degree harmonics from OSU89B were also considered. A 
number of numerical experiments were performed that lead to the following conclusions: 


(a) The currently available global geopotential model OSU89B alone is expected to yield 
geopotential differences between stations separated by 30°, accurate to about 86 kgal cm. 
The future models (augmented by OSU89B) can improve this accuracy to about 
81 kgal cm (TOPEX/OSU 89B) and 74 kgal cm (GPB/OSU89B) respectively. 

(b) Introduction of gravity disturbance measurements in terrestrial 2° caps 
reduces the previous error estimates to the following: 23 kgal cm (OSU89B), 17 kgal cm 
(TOPEX/OSU 89B) and 12 kgal cm (GPB/OSU89B), when pessimistic error estimates are 
used for the gravity disturbance measurements (mo — 2 mgal). With mo = 1 mgal the case 
GPB/OSU89B yields an error of about 9 kgal cm for 30° station separation. The error 
estimates for these cases were computed using both truncation theory and lsc (ring 
averages) and the results from the two techniques were compared. It was found that the 
lsc error estimates are always smaller than the ones obtained from truncation theory, as 
mandated from theory. The largest difference between the two error estimates was found 
to be about 21%. 

(c) When gravity disturbance data at the altitude of GP-B (about 600 km) were 
introduced, a moderate (7%) improvement in accuracy, over the corresponding case 
without such data, was found. In both cases, the reference geopotential model used was 
complete to degree 45, obtained from the analysis of the GPS tracking data on GP-B. 
However, gravity disturbance data at this altitude are unable to resolve medium and high 
frequency variations of the gravity field and thus the result in this case is inferior by about 
5 kgal cm to the result obtained from the combined GPB/OSU89B high-degree model 
(complete to degree 250). 

To enrich the data at altitude with more high-frequency information, it is 
recommended here that additional measurements of a higher-order gradient of the 
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disturbing potential made from a lower flying spacecraft, be incorporated in the estimation. 
In this direction, the gradiometer data from ARISTOTELES mission can provide a 
significant contribution. In addition, it should be emphasized that the error estimates 
reported here correspond to a "worst-case" scenario where only one pair of benchmarks is 
considered for the intercontinental connection. Additional benchmarks on each continent, 
connected with leveling lines, can provide a better network configuration and yield an 
improvement on the accuracy of the intercontinental connections up to 25%, as the study 
by Hajela (1983) has indicated. 

Finally, the results reported here are promising enough to warrant an actual testing 
of the technique. For this purpose, stations whose geocentric coordinates are accurately 
known (e.g. SLR sites or VLBI stations connected to a geocentric system using GPS) and 
between which the geopotential difference has been estimated independently using spirit 
leveling and gravimetry can be used as test sites. At present (1991), collection of gravity 
disturbance measurements (using relative GPS positioning and gravimetry) in caps 
surrounding these test sites, will enable testing of the procedure described in section 4.2. 
As it was discussed in that section, a cap size of 2° and an approximate spacing of 6' 
between the points where the gravity disturbances are determined, are optimum parameters 
for the observational setup, provided a state-of-the-art high-degree global geopotential 
model (e.g. OSU89B) is used as reference. In the actual implementation, least-squares 
collocation using the original gravity disturbance data (as opposed to ring-averages) 
should be used to maintain highest computational rigor. In addition, detailed elevation 
information around the test sites should be used for the consideration of the terrain effects 
by means of analytical continuation, as it was discussed in section 2.3. 

In a more future time frame (1995), the availability of the data from the anticipated 
satellite missions (TOPEX/Poseidon, Gravity Probe-B, ARISTOTELES), will enable to 
improve the above scheme in two ways. First, by the use of global geopotential models 
with more accurately estimated lower-degree harmonics than those of OSU89B, as such 
models will become available from the analysis of the global sets of observations collected 
by these missions, and second by the use of gravitational information in caps at altitude 
(GP-B, ARISTOTELES), as discussed in section 4.3. 



APPENDIX A 


RECURRENCE RELATIONS FOR THE TRUNCATION COEFFICIENTS IMPLIED 
BY PIZZETTI’S EXTENSION OF HOTINE’S KERNEL 


The development of recurrence relations for the coefficients Vo), defined by: 

Qi (H, Vo) = f H(R/ r , v)P n (cosv) sinvdv ; r ^ R (A.l) 

'Vo 

is based on recurrence relations of the Legendre polynomials. Denoting P n (t) (or simply 
P„) the n* - degree Legendre polynomial, and Pn(t) (or simply Pn) its first derivative 
with respect to the argument t, the following relations hold true (Hobson, 1965, Ch. II, 


Sect. 20): 

nP„ = -(n-l)P n -2 + (2n-l)tP„-i (A.2) 

(2n + 1)P„ = P n +i - Pn-1 (A.3) 

(l-t2)Pn = n(P n -l-tPn) (A.4) 

Pn(" 1) = (- D n ( A - 5 > 

The following notation is introduced: 

k = y ; 0 <k £ 1 (A.6) 

t = cosv ; -l£t<l (A.7) 

D 2 = 1 - 2kt + k 2 (A.8) 


Accordingly, Pizzetti's extension of Hotine’s kernel takes the closed form: 
H(k,t) = ^-lnp±J^) ;t*l 
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and the series expansion form: 


oo 



(A. 10) 


With the change of integration variable from \|/, to t = cosy, and the substitution 
x = cos\j/ 0 , equation (A.l) becomes: 


-£ 


Qn(k 

or, due to (A.9): 


H(k,t) P n (t)dt 


Qn(k,x) = 2k j 5^dt - j ln(D±jL^)p n(t)dt 


(A.l 1) 


(A. 12) 


It will be understood in all subsequent derivations that \|r 0 * 0 so that x * 1. The case 
\j/o= 0 will be considered separately at the end. 


For the purposes of the subsequent derivations, a number of auxiliary quantities 
are defined next: 


Into = 



Pn(t)dt 


(a) 
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i;(k,x) = 


K>,x) = 


M n (k,x) = 






(b) 




(A. 13) 
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To develop recurrence relations for Qn(k, x), one starts with recurrence relations for 
some of the auxiliary quantities defined before. 

From (A.2) one has: 

Pn(x) = l [(2n - 1) xP n .i(x) - (n-l)P n - 2 (x)] ; n ^ 2 (A.14a) 

which, with starting values: 

P 0 (x) = 1 ; Pi(x) = x (A. 14b) 


establishes a recurrence relation for the Legendre polynomials. 

From (A.3) and (A.5) one has: 

In(x) = 2nVr [Pn + 1 (x) - Pn ’ l(x)] ; n ^ 1 (A15) 

Although equation (A. 15), along with the recurrence (A.14a,b), are sufficient to evaluate 
I„(x), an independent recurrence for I n (x) may be developed as follows. Equation (A. 15) 
is written as: 


- 01 - 2)I„ . 200 = - £^P» . l(x) + . 3 (x) 

In addition, from (A. 15) one has: 

(2n - l)xl n . i(x) = xP„(x) - xP n . 2 (x) 

Now, from equation (A.2) one may write: 

xP»(x) = JL±l- p, , tl « +s ^ T l > n ., ( x) 


(A. 16) 


(A. 17) 


(A. 18a) 


and: 

- xPn - 2 (x) = - 2 ^' ^Pn - l(x) - 2 ^T^P n - 3 W 


(A. 18b) 
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Adding by parts (A.16), (A. 18a) and (A.18b), and taking into account (A.17), one has: 

(2n - l)xl n . i(x) - (n - 2)I n . 2 (x) = n + \ [P n + i(x) - P n - i(x)] 

zn + 1 

which, due to (A.15), finally establishes the following recurrence for I n (x): 

I n(x)=-h-[( 2n -l)x I „.i(x)-(n-2)I n . 2 (x)] ; n>2 (A. 19a) 

n + 1 

with starting values: 

I 0 (x)=l+x ; Ii(x)=i(x2- 1) . (A. 19b) 

For Lil(k, x), due to (A. 3), one has: 


li(k, x) = 



dt 


so, due to (A. 13d): 

L' n (k, x) = -- - l - [M n + i(k, x) - M n -i(k, x)] 

zn + 1 (A.20) 

With dv = Pn(t)dt (v = P n (t)), and w = D' 1 (dw = kD‘ 3 dt), integration by parts yields for 
M„(k, x): 

M n (k,x)=^ - ^k) - kK n(k, x) (A.21) 


where: 

D x = (1 - 2kx + k 2 ) Vz 


(A.22) 


Equations (A.20), (A.21) and (A.15), yield for L„(k, x): 


Consider K„(k, x) next. The definition (A. 13c) and equation (A.2) imply: 
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(n + l)< + 1 +nK n . 1 =(2n + l)J^-^dt 
However, it can be easily verified that: 


t _ 1 + k 2 T 1 1 

D 3 2k [d 3 D(l+k 2 ) 


(A.24) 


(A.25) 


Due to (A.25), the right-hand side of (A.24), becomes: 

//i . I \ |* tP n (t), ^ , i\l + k 2 ^ 2n + 1 ( ^*n(f 
(2n+l)J i -^ r <it-(2n+l)-^ jr -K n - ^ J D 


so that, taking also into account the definition of L^(k, x), equation (A.24) becomes: 


2k[(n + 1)K^ + j + nK^ . J = (2n + 1)[(1 + k 2 )K n - lJ 


(A.26) 


Eliminating from equations (A.23) and (A.26), one obtains: 
kK^ + : (k, x) - (1 + k 2 )K n (k, x) + kK^. j(k, x) = - ^ 


In addition, direct integration yields: 


fiiL = _l_ 
J D3 kD 


; n £ 1 (A.27) 


-idL = — L(1 - kt + k 2 ) 
D 3 k 2 D 


Hence, a recurrence relation for K n (k, x) may be established as follows: 


(A. 28) 


;<k,x) = l^i<,. 1 (k,x)-K„. 2 (k,x)-a j ^ ; n - 2 


In-l(x) 



with starting values: 
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K 0 (k, x) = 


kD 


K'(k, x) = 



kt + k 2 )| 


'-l 


(A. 29b) 


One may now express Kn+i(k, x) in terms of K n (k, x) and K n -i(k, x) from equation 
(A.29a), and then substitute this expression in equation (A.23), so that finally the 
following recurrence relation is established for L^k, x): 


L,i(k, x) = 


1 

2n+ 1 


I„(x) 


2kK n l (k, x) - (1 + k 2 )K;(k, x) + 2(n + 1>^ 


; n > 2 (A. 30a) 


with starting values: 


L 0 (k, x) = - C. 


X 


-1 


(k, x) = -J2-(l+kt 
3k 2 



(A.30b) 


which were obtained by direct integration applied to the definition of L^(k, x) given in 
(A. 13b). Notice that the recurrence (A.30a) may actually start from n = 1, however, it is 
computationally convenient to start it from n = 2, the same degree from which the P n (x), 
I„(x), Kn(x) recurrence relations start. 

With the previous derivations as a background, one proceeds now with the 
development of recurrence relations for Qn(k, x), as follows. From equations (A. 12) and 
(A. 13b) one has: 

Qn(k, x) = 2kL>, x) - j In ( P + fc-t) P n (t)dt . 

Define: 

R'„(k, x) = - j ln(Q±^)p n (t)dt (A.31) 


so that the previous equation becomes: 
Qn(k, x) =2kl4(k, x) + R^(k, x) 


; n ^0 


(A. 32) 



Due to equation (A.2), equation (A.31) becomes: 
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r' - n - 1 p ' 2n - 

K n~ n K n-2 


which due to (A.4) becomes: 






(A.33) 


The integral on the right-hand side of (A.33) may be evaluated using integration by parts. 
Let: 


dv = Pn-i(t)dt => v = P n -i(0 


and: 


w = (l-t 2 )ln( D [ k t - t ) 

It is a matter of algebraic manipulations, to show that: 
dw = [- 2tln( P [jy 1 ) + d ’ d + l ] dt 

Hence, integration by parts in equation (A.33) yields: 


where: 


a>Jv p "( t)dt 


(A.34) 


(A. 35a) 


(A.35b) 


(A.36) 


(A.37) 


and: 


B 'n = j ^^7^) tP n 


(t)dt 


(A.38) 
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For B n , due to equation (A. 2) and the definition (A.31) one has: 




For An, due to the definition (A. 13b) one has: 

A '„-l = k j ^P„-l(t)dt-L„_,, 
which due to equation (A.2) becomes: 


(A.39) 


An-1 = N 


n ( p n (t) d n- 1 f Pn-2 (t) 
2n- 1 J D 2n- ljj D 


dt 


"^nt 


which finally becomes: 


2n - 1 / J—L— x* +1 l J + 2n- 1 l , 

n(n - 1) n l Ln - 1 ^ n ^- 2 \ n(n . 1} n-i 


(A.40) 


Substituting (A.39) and (A.40) in equation (A.36), one obtains the following recurrence 
relation for Rn(k, x): 

R'„(k, X) = j(n - 2)(n - 1) R„_ 2 (k, x) + (2n - 1)(1 - x 2 )lnj^i^j P„.,(x) 


-k 


+ 


nL n (k, x) + (n - l)L n2 (k, x) ] 
(2n-l)[L n . 1 (k, x)-I n .i n > 2 


(A.41) 


which requires starting values for R^(k, x) and Ri(k, x). These may be obtained by 
direct integration applied to equation (A.31) that defines Rn(k, x). Let: 



r i(k, t)= -Jtln p^-^ dt 
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It can be shown that: 


r^k, t ) = -Lib 2 ln p - tJL zl) - (1 + k) 2 ln(D + k + 1) + (1 - k) 2 ln(D - k + 1) + 2kD 
21c 1 — t ' 


(A.42a) 


r'(k, t) = ^{k[ln(D+k+l)-ln(D-k+l)]-|[l|--D 2 + (l-k 2 )D] 

-kt 2 ln P|k-* )) 


(A.42b) 


so that the starting values R^(k, x) and R{(k, x) are given by: 


R 0 (k, x) = r 0 (k, t) | t 


Rj(k, x) = rj(k, t) 


(a) 


(b) 


(A.43) 


The recurrence relation for Rn(k, x) was the last relation required to establish a 
recursive algorithm for the computation of the truncation coefficients Qn(k, x). In 
summary, for y 0 * 0 (x = cos\jr 0 * 1), the recursive evaluation of Qn(k, x) is 
accomplished by computing the following quantities: 

1 . P n (x) from equations (A. 14a, b) 

2 . I n (x) from equations (A. 19a, b) 

3 . Kn(k, x) from equations (A.29a, b) 

4. Ln(k, x) from equations (A.30a, b) 

5 . Rn(k, x) from equations (A.41) through (A.43) 

6. Finally Qn(k, x) is given by equation (A.32) 
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Two issues concerning the above recurrence should be mentioned: 

a) The recurrence for Kn(k, x) is only required for the evaluation of L^(k, x). In the 
special case that r = R (i.e. k = 1), Jekeli (1979) has succeeded in eliminating Kn(k, x) 
from the recursive algorithm, by manipulating equations corresponding to (A.29a) and 
(A.30a). However, this technique cannot be applied in the current case where k assumes 
arbitrary values in the interval (0, 1]. 

b) Although mathematically correct, and useful for low altitude applications, the 
recurrence presented here for Qn(k, x) has numerical instability problems for high altitude 
cases. The instability is introduced to the algorithm through the use of equation (A.29a) 
that determines Kn(k, x). As it can be seen from (A.29a), the division by k introduces 
numerical problems as k->0 (i.e. as the altitude increases). Through numerical tests 
(using double precision arithmetic) it was verified that for altitude equal to 600 km, and 
\|r 0 as small as one degree, the coefficients Qn(k, x) become unreliable after n = 600. 
This problem is of the same nature with the numerical instability encountered by 
Shepperd (1979) in the recurrence relation which he derived for the altitude generalized 
truncation coefficients corresponding to Stokes' kernel (ibid, p. B-3). 


Finally , to complete the determination of Q n (k, x), one needs to derive the 
expression for Qn(k, 1), i.e. consider the case \jr 0 = 0 (x = 1), which was excluded from 
the above algorithm. Although H(k, t) has a singularity at t = 1 (see equation (A.9)), 
Qn(k, x) is nevertheless well-defined for x = 1. Taking into account equation (A. 10), it 
can be easily seen that: 


<Wk. 1 ) = ^r k ” + 1 


n > 0 


(A. 44) 



APPENDIX B 


TRUNCATION COEFFICIENTS IMPLIED BY THE 
KERNEL H*(k, t) = H(k, t) - k - 1 k 2 t 


In case the gravity disturbance, 5g(r, 0, X), contains no zeroeth- and first-degree 
terms, one needs to determine the truncation coefficients Qn (k, x) implied by the kernel 
function H*(k, t), where: 

H*(k,t) = H(k,t)-k-|k 2 t (bj) 

or, in series expansion form: 


H*(k, t) = X 2n TT- k " +lp n( t ) • 

n=2 n+1 

The truncation coefficients Qn (k, x) are defined by: 
Q^(k,x)=| H*(k,t)P n (t)dt 

or, 

Qn(k, x) = Qn(k,x)-k| P n (t)dt-|k 2 | tP„(t)dt 


(B.2) 


(B.3) 


(B.4) 


Utilizing the results of Appendix A, one has from (A. 13a): 


Pn(t)dt = I n (x) 

1 (B.5a) 
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while, due to equation (A.2): 

/ i t, ( .^^rv lW+ -T , »., W 

The last equation, due to (A. 19a), may be rewritten as: 

Due to (B.5a,b) equation (B.4) finally becomes: 

Qn(k, x) = Qn(k, x) - jkl n (x) + 2 ^ 4 [(n + l)xl n (x)+ I„.i(x)] j ; n > 2 (B.6) 


and, along with the starting values: 
Qo(k, x) = Q 0 (k, x) - 
Qi(k, x) = Qi(k, x) - 


k(x + 1) + (x 2 - l)j 

(a)| 

| (x 2 - 1) + ^ (x 3 + l)] 

(b) / 


(B.5b) 


(B.7) 


(which have been obtained by direct integration performed in equation (B.4)) defines the 
recurrence relation required to evaluate the coefficients Qn(k, x), from the already derived 
recurrences for Qn(k, x) and I n (x) (see Appendix A). 

The recurrence relation (B.6) is valid regardless of the value of \|/ 0 . If \/ 0 = 0, 
then Qn(k,l) in (B.6) are obtained from equation (A.44); otherwise the Qn(k, x) are 
obtained from equation (A.32). If one rewrites equation (B.4) as: 

(X(k, X) = Qn(k, X) ■ kj Po(t)P„(t)dt -jk 2 j P,(t)P„(t)dt (B.8) 


which is valid for n > 0, then it can be seen easily that: 

(a) (£(k, 1) = Qn(k, 1) for every n > 2 

(b) Qo(k, 1) = Q*(k, 1) = 0 , while Q£(k, x) and Q*(k, x) are not equal to zero for x* 1. 



APPENDIX C 


DEGREE VARIANCES OF A FIRST-ORDER GAUSS-MARKOV 
PROCESS ON THE SPHERE 


Consider the homogeneous and isotropic first-order Gauss-Markov stochastic 
process on the sphere, defined through its covariance function: 

o(\|/) = c.e x v ; c>0,X>0 (C.l) 


where \j f is the spherical distance in radians. The degree variances of the above 
covariance function constitute the power spectrum of the process, and are given by the 
Legendre transform of o(y). Using unnormalized spherical harmonics, due to isotropy, 
one has for the degree variance at degree n: 

J r 71 / 2 k 

j a(Y)P n (cos\|r)sin\|/d\j/da 

V=0 Ja=0 (C.2) 

or, due to (C.l): 

P 

o n = 2n + - c I e‘^vp n (cos\j/)sin\jtdY 

2 Jo (C.3) 


Let: 


ajM= exp(->.y)P n (cosY)sin\j/dY 

Jo (C.4) 


so that c n is given by: 

o n = 2n + -LcA n ()i) 


(C.5) 
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One proceeds next with the derivation of an expression for A n (X). Changing the variable 
from \|/ to t = cos\|/, one has from (C.4): 


a„M = exp(-Xcos' 1 t)P n (t)dt 
which, due to equation (A. 3) of Appendix A, yields: 

expC-teos-UXP^j - P n l )dt 


AnW W! 


(C.6) 


With w = exp(-Xcos _1 t) (hence dw = Xw (1 - t 2 )' I/2 dt), and dv = (P n '+i - P n '-i)dt (hence 
v = Pn+l(t) - Pn-l(t)), integration by parts in equation (C.6), yields: 


A nW = 2^Vl| exp(-Xcos' 1 t)(l - 1 2 Y k (P n .i - P„ + i)dt 


(C.7) 


where, in deriving (C.7), use was made of equation (A.5) of Appendix A. From 
Hobson (1965, equation 20.36), one has: 


(2n + l)(t 2 - l)P n = n(n + l)(P n+1 - P n .j) 
according to which, equation (C.7) becomes: 

A„W = exp(-Xcosh)(l - t 2 fhp n dt 


(C.8) 


(C.9) 


Now let: 


w 


= exp(-Xcos 4 t)(l - 1 2 ) 1 ^ 2 


so that: 


dw = exp(-Xcos _1 t)[x - t(l - t 2 )^ 2 ]dt 


and, dv = P„ dt (hence, v = P n (t)). Then, integration by parts in equation (C.9) yields: 
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, ... f expt-XcoS'Ujtl - t^'^tPndt - XA n (x) 

n ( n + !)[Ji 

which, due to equation (A.4) of Appendix A, yields: 

AnW = n (n V l) expC-Xcos-U^l - 1 2 \ ^ 2 P n -idt 

- -ij expC-teos-^l - t 2 )^ 2 P n dt - XAj(x) 
which, taking into account equation (C.9), reduces finally to: 

[(n + if + >. 2 ]A n (>.) = xj expC-Xcos 1 ^! - t 2 )^P n .idt 


(C.10) 


Evaluating (C.10) for n = n - 2 and subtracting (C.10) from the corresponding equation, 
one obtains: 

[(n - if + k 2 ]A n . 2 M - [(n + if + A. 2 ]A n (x) * xj expC-Acos-^l - t 2 ) - ^(P n -3 - P n -i) dt 

The right-hand side of the last equation equals, due to (C.7), to (2n - 3) A n _ 2 (X) so that 
the preceding equation, after the algebraic simplifications, yields: 



(C.ll) 


Considering now equation (C.5), it can be easily seen that the degree variances of a(\|/) 
may be evaluated from the following recurrence relation: 


_ 2n + 1 
2n - 


2 + (n - 2f 
2 + (n + 1) 


n £ 2 


(C.12) 


which requires starting values for o Q and cti. These are obtained by direct integration 
performed in equation (C.3) which yields: 
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(a) 

2 r + 1 


(C.13) 
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